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The paper gives full details of the computation within the canonical formalism of Arnowitt, Deser, 
and Misner of the local-in-time part of the fourth post-Newtonian, i.e. of power eight in one over 
speed of light, conservative Hamiltonian of spinless compact binary systems. The Hamiltonian 
depends only on the bodies’ positions and momenta. Dirac delta distributions are taken as source 
functions. Their full control is furnished by dimensional continuation, by means of which the 
occurring ultraviolet (UV) divergences are uniquely regularized. The applied near-zone expansion of 
the time-symmetric Green function leads to infrared (IR) divergences. Their analytic regularization 
results in one single ambiguity parameter. Unique fixation of it was successfully performed in 
T. Damour, P. Jaranowski, and G. Schafer, Phys. Rev. D 89, 064058 (2014) through far-zone 
matching. Technically as well as conceptually (backscatter binding energy), the level of the Lamb 
shift in quantum electrodynamics is reached. In a first run a computation of all terms is performed in 
three-dimensional space using analytic Riesz-Hadamard regularization techniques. Then divergences 
are treated locally (i.e., around particles’ positions for UV and in the vicinity of spatial infinity for IR 
divergences) by means of combined dimensional and analytic regularization. Various evolved analytic 
expressions are presented for the first time. The breakdown of the Leibniz rule for distributional 
derivatives is addressed as well as the in general nondistributive law when regularizing value of 
products of functions evaluated at their singular point. 


I. INTRODUCTION 

In the early days of Einstein’s theory of gravity, the 
theory of general relativity as Einstein coined it, the 
weak-field and slow-motion expansion of the field equa¬ 
tions, nowadays known as post-Newtonian (PN) expan¬ 
sion, played a crucial role already. The first post- 
Newtonian (1PN) approximation revealed a convincing 
understanding of the perihelion advance of Mercury [H . 
Later on, in 1919, the light bending at the limb of the 
Sun was measured in agreement with the 1PN predic¬ 
tion of Einstein’s theory i 1- A rich 1PN scenario 
came up with the discovery of the Hulse-Taylor binary 
pulsar PSR B1913+16 in 1974, where after four years 
of observation even the dissipative two-and-a-half post- 
Newtonian (2.5PN) effect could be measured in agree¬ 
ment with general relativity in the form of energy loss 
through gravitational radiation damping [4]. Also the 
precession of the proper angular momentum (spin) of the 
pulsar PSR B1913+16 from spin-orbit coupling, a 1PN 
or 1.5PN effect according to the counting of the spin as 
of 0PN or 0.5PN order, respectively, could be seen 
The discovery of the double pulsar system J0737—3039 
in 2003 allowed measurements of even more 1PN effects 
than with the Hulse-Taylor binary pulsar j(| . Also in the 
solar system and in the gravitational field of the Earth 
the 1PN approximation is crucial [ 3 , 0 |- In the observa- 
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tion of the Hulse-Taylor pulsar, even 2PN measurement 
accuracy has been achieved in the periastron advance [9j], 
deduced for the first time in M with full details given in 
0. In the long-term run, the observation of the double 
pulsar system is likely to allow detailed measurements of 
conservative 2PN and dissipative 3.5PN effects ftii. IT3 - fl4| . 
On the other side, in the coming years gravitational-wave 
astronomy will come to operation through advanced (or 
second-generation) LIGO and Virgo observatories fl5Ulfil| 
and the cryogenic KAGRA detector (formerly known as 
LCGT) [13| . Even more sensitive third-generation de¬ 
tectors, like ET [l8|, are already under consideration. 
Then still higher-order PN effects will become important. 
Great efforts have already been undertaken in the precise 
data analysis of possible wave forms (see, e.g., fl9H22| ). 
and as well as in the computation of the orbital and spin 
dynamics through higher PN orders including the gravi¬ 
tational wave emission (see, e.g., [23l - l30i ]). 

The complete conservative third post-Newtonian 
(3PN) binary dynamics has been fully achieved for the 
first time in 2001 [Hj, based on earlier work [32l - [35[ , 
rooted in the canonical formalism of Arnowitt, Deser, 
and Misner (ADM) with application of a coordinate sys¬ 
tem of maximal isotropy in three-dimensional space 36]. 
Around 2004 and later in 2011 the final results of [3ll — 
[3|| became fully confirmed pul - Eol] . in using completely 
different techniques but yet all with the employment of 
harmonic coordinate conditions. Let us stress that the 
rather tedious (but much easier than reported in the 
present paper) calculations performed in pilUdbl ] were 
not just a matter of straightforward computational ef- 
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forts. The standard regularization Hadamard “partie- 
finie” technique together with the Hadamard “partie- 
finie” based Riesz analytic continuation supplemented by 
the Schwartz distribution theory in (3+l)-dimensional 
space-time resulted in ambiguities [32|, [ 33 , which yet 
could be parametrized by two variables only, w s and Wk 
[34 ] (for a detailed presentation of the mentioned tech¬ 
niques, see jU| and Appendix El of the present paper). 
Matching to the Brill-Lindquist initial value solution for 
a two-black-hole system [421 ], a solution which was shown 
derivable from “fictitious” point-mass sources of Dirac 
5-function type (one for each black hole) ]4J, brought 
out (by the very definition of w s ) ui s = 0 [33, [HJ, and 
implementation of Lorentz invariance in the form of ful¬ 
fillment of the Poincare algebra resulted in Wk = 41/24 
(see [35] and Sec. IX of the present paper). Another pro¬ 
cedure to obtain w s = 0 could have been to insist on the 
“tweedling of products” structure in the course of com¬ 
putation of re gula rized value of the product of singular 
functions (see [33|, [44[ and Appendix IA4I of the current 
paper) or, connected herewith, by taking just the finite 
terms (without limiting procedure) in the two-body re¬ 
striction of the many-body static potential derived in 
[45[ . Only by making use of the technique of dimen¬ 
sional regularization a uniform treatment was achieved 
in Ref. [3l|, with confirmation of the numerical values 
for w s and Wk found earlier (a summary of the applied 
dimensional-regularization techniques is presented in [40] 
and in Appendix |A] below). For extensive application of 
the dimensional regularization in the generalized ADM 
formalism for spinning compact binaries the reader may 
consult [28j . 

It might be worth comparing the derivation of the 
3PN ADM Hamiltonian in [31H35I ] with the history of 
the harmonic-coordinate-based calculations, which led 
to the results achieved in [39], where finally also di¬ 
mensional regularization played the crucial role. In 
the papers fil], ®] a manifest Lorentz-invariant “ex¬ 
tended” Hadamard regularization procedure was devel¬ 
oped, which allowed the calculation of the 3PN binary 
dynamics with only one ambiguity parameter A. The 
comparison, for circular motion, of energy in terms of or¬ 
bital angular frequency (which is a coordinate-invariant 
or gauge-invariant relation) revealed the relation A = 
—3w s /ll — 1987/3080 Eflt. Evidently, the application of 
the Poincare algebra in [35| is the equivalent of the man¬ 
ifest Lorentz invariance in [id], HU- However, as shown 
in [39|, the extended Hadamard regularization proce¬ 
dure, being incompatible with Schwartz distribution the¬ 
ory [see, e.g., Eqs. (7.17) and (8.34) in [13], which should 
be compared with Eqs. (IA90bl) and (IA92bl) . respectively, 
of our paper], could not be made compatible with di¬ 
mensional regularization [see, e.g., Sec. Ill D in [39j|]. 
Therefore all terms of the extended Hadamard procedure 
different from the standard Hadamard (“pure Hadamard- 
Schwartz” in HI) ones, had to be traced back to their 
origins and eliminated [see, e.g., Eq. (3.55) in [39|] before 
dimensional generalization and regularization could be 


employed. There exists another derivation of the 3PN bi¬ 
nary dynamics in harmonic coordinates made in [4fi| , us¬ 
ing an effective field theory approach with space-Fourier 
transformed fields, which was from the very beginning 
put on a (d+l)-dimensional space-time footing with di¬ 
mensional regularization. 

It is interesting to note that the still other derivation 
of the 3PN binary dynamics in [37, HU, based on the 
Einstein-Infcld-Hoffmann surface integral technique (50| , 
is a purely (3+l)-dimensional one. The used surface- 
integral technique allowed a manifest Lorentz-invariant 
calculation with surface integrals defined in the smooth 
vacuum regime. Divergent integrals entered only on tech¬ 
nical reasons in simplifying computations. (In the paper 
by Einstein, Infeld, and Hoffmann [50] the field singu¬ 
larities were already seen in correspondence with Dirac 
5-functions, though their treatment as full field sources 
had had to await for later developments in the form of 
Infeld’s “good 5-functions,” see, e.g., the book by Infeld 
and Plebanski [44j.) Remarkably, the Brill-Lindquist ini¬ 
tial value solution in (3+l)-dimensional space-time men¬ 
tioned above has originally been obtained from the pure 
vacuum Einstein field equations without facing any phys¬ 
ical or geometrical divergences [42j. 

The fourth post-Newtonian (4PN) conservative dy¬ 
namics of two-point-mass systems has been completed 
only quite recently 5l|, based on previous calculations 
[53 . 153 ]. The results of Refs. [5jjj53| were in part con¬ 
firmed by [53 - 1571 (see also 0,HU)- The present paper, 
announced in [51[, delivers the details of the involved 
computations of the 4PN conservative two-point-mass 
ADM Hamiltonian, where ultraviolet (UV) and infrared 
(IR) divergences have to be tackled simultaneously and 
where a mixed dimensional-analytic regularization treat¬ 
ment is needed to be applied (like in [60|). Pure dimen¬ 
sional regularization allows one to control UV divergences 
only and Refs. [52;. _53| have uniquely regularized all UV 
divergences. Regularization of IR divergences turned out 
to be ambiguous (what is discussed in the present paper) 
and this ambiguity was resolved in Ref. 151| by taking 
into account the breakdown of a usual PN scheme (based 
on a formal near-zone expansion) due to infinite-range 
tail-transported temporal correlations found in Ref. [6JL|. 
Reference [5l| showed that the total 4PN conservative 
Hamiltonian is the sum of instantaneous (local-in-time) 
near-zone Hamiltonian and time-symmetric but nonlocal- 
in-time tail Hamiltonian. The present paper is devoted 
to details of computation of the 4PN near-zone local-in¬ 
time Hamiltonian. A rather nontrivial application of Ref. 
[51 ] to 4PN-accurate generalization of effective-one-body 
approach to compact binary dynamics has been worked 
out in Ref. already. 

Since the seminal work by’t Hooft and Veltnran [63j . 
dimensional regularization (for history see H) has be¬ 
come very popular in quantum field theory (see, e.g., 
[65|'l. Even the famous Lamb shift with simultane¬ 
ously occurring UV and IR divergences found its ele gan t 
computation by applying dimensional regularization [66(. 
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However, whereas UV divergences are nicely controlled 
by dimensional regularization, IR divergences can pose 
problems (see, e.g., Refs. [671,[ 63 ). This also happens in 
our approach with the final solution indicated above. It 
is worth pointing out that also the Lamb shift calcula¬ 
tion of Ref. [66] shows up an undefined constant in the IR 
sector, which gets fixed by some dimensional matching. 

A short overview of the various sections will give the 
reader some help in orientation. Section II presents the 
ADM formalism for two-point-mass systems developed in 
(d + l)-dimensional space-time. It also gives the tran¬ 
sition to the Routhian functional which is crucial for 
the obtention of Hamiltonian which depends on matter 
variables only. Section III expands the nonpropagating 
part of the ADM structure for two-point-mass systems 
through 4PN order. It ends with the formulas for the 
4PN-accurate reduced Hamiltonian. Section IV is con¬ 
cerned with the derivation of field equations for the prop¬ 
agating degrees of freedom valid to the next-to-leading 
order (which is enough for our purposes). It also in¬ 
cludes the formal near-zone (and thus PN) expansion of 
the time-symmetric solution of field equations. Section 
V develops the 4PN-accurate Routhian functional, with¬ 
out expanding propagating degrees of freedom into the 
PN series. Section VI presents the 4PN-accurate conser¬ 
vative Hamiltonian dependent on matter variables only, 
with still non-PN-expanded propagating degrees of free¬ 
dom. Section VII gives the 4PN-exact near-zone conser¬ 
vative matter Hamiltonian with fully PN-expanded prop¬ 
agating degrees of freedom. Section VIII delivers some 
details of the computation and regularization of integrals, 
taking into account the various appendices, and presents 
the fully explicit result for the total 4PN-accurate mat¬ 
ter conservative Hamiltonian. Section IX is devoted to 
checks with the aid of the Poincare algebra. All technical 
calculational details are shifted to three appendices. Ap¬ 
pendix A is devoted to all regularization procedures used. 
Appendix B gives a variety of needed inverse Laplacians. 
Finally, Appendix C presents field functions in fully ex¬ 
plicit forms. 

We employ the following notation: x = (x l ) (i = 
1,... ,d) denotes a point in the d-dimensional Euclidean 
space R d endowed with a standard Euclidean metric and 
a scalar product (denoted by a dot). Spatial latin indices 
run through 1 ,... ,d, and space-time greek indices vary 
from 0 to d. Letters a and b (a, b = 1,2) are particle 
labels, so x a = (x®) £ R d denotes the position of the 
ath point mass. We also define r a := x — x a , r a := |r a |, 
n a := r a /r a ; and for a ^ b, r ab := x a - x b , r a b ■= |r a6 |, 
n ab '■= r ab/i'ab'i | • | stands here for the Euclidean length of 
a vector. The linear momentum vector of the ath particle 
is denoted by p a = (p a i), and m a denotes its mass pa¬ 
rameter. An overdot, as in x a , means differentiation with 
respect to time coordinate t. The partial differentiation 
with respect to x 1 is denoted by di or by a comma, i.e., 
difj) = and the partial differentiation with respect to 
x® is denoted by d a i- We abbreviate Dirac delta distri¬ 
bution <5(x — x a ) by d a (both in d and in 3 dimensions); 


it fulfills the condition f d d x5 a = 1. Flat d-dimensional 
Laplacian is denoted by A^, whereas A (without any 
subscript) is reserved for Laplacian in d = 3 dimensions. 
Extensive use has been made of the computer algebra 
system Mathematica. 

II. THE ADM CANONICAL FORMALISM 
IN d SPACE DIMENSIONS 

We consider here a system of two point masses, i.e. 
monopolar, pointlike bodies, which interact gravitation¬ 
ally according to general relativity theory. We model 
point masses by means of Dirac delta distributions S a - 
Let D := d + 1 denote the (analytically continued) space- 
time dimension. The ADM canonical approach [36[ uses 
ad + 1 split of the coupled gravity-plus-matter dynam¬ 
ics and works with the pairs of the canonical variables: 
positions x a and momenta p a of point masses, and for 
the gravitational field the space metric 7 ^ = < 7 ^ induced 
by the full space-time metric gp v on the hypersurface 
t = const, and its conjugate momentum 7 t ij = tt ji . 

The full Einstein field equations in D dimensions in an 
asymptotically flat space-time and in an asymptotically 
Minkowskian coordinate system, written in the canonical 
variables introduced above, are derivable from the Hamil¬ 
tonian ( 6 fj (in units where 167 t Gd = c = 1 , with Gd de¬ 
noting the generalized Newtonian gravitational constant 
and c the speed of light) 

H = J d d x (. NU-N l Ut) 

+ <f d d ~ 1 S i dj{~/ij - Sijjkk), ( 2 . 1 ) 

Ji° 

where N and N' 1 are called lapse and shift functions and 
the super-Hamiltonian H and super-momentum Hi are 
defined as follows: 

H(x a , p a , 'Ytj, n ij ) := ^N 2 (T 00 - 2G 00 ) , (2.2a) 

n t (x a , p a , := fyN (T° - 2G°) . (2.2b) 

Here T^ v and G^ denote the energy-momentum and the 
Einstein tensor, respectively, and 7 is the determinant of 
the d-dimensional matrix ( 7 ^). In the surface integral 
in m *° denotes spacelike infinity and d d 1 5 1 ,; is the 
(d — l)-dimensional out-pointing surface element there. 
In terms of lapse and shift functions, the space-time line 
element reads (we employ signature d — 1 , x° = t) 

ds 2 = g^ v dx^dx v 

= -N 2 {dx 0 ) 2 + 'y ij (dx i + N i dx°)(dx j + N j dx°). 

(2.3) 

The lapse and shift functions are Lagrangian multipli¬ 
ers only and deliver the famous Hamiltonian and momen¬ 
tum constraint equations of the Einstein theory, 

H = 0 , Hi = 0 . 


(2.4) 
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The fulfillment of the constraint equations together with 
d+1 additional (independent) coordinate conditions, see 
Eqs. (El below, reduces the total Hamiltonian of Eq. 
EH) to the reduced one, H le d, which obviously is a pure 
surface expression. This is the ADM Hamiltonian which 
evolves the independent dynamical degrees of freedom of 
the total system. Here the lapse and shift functions are 
not involved at all, only the space-asymptotic numerical 
value 1 of the lapse function enters. The lapse and shift 
functions are to be calculated by making use of the coor¬ 
dinate conditions within the remaining set of the Einstein 
field equations. 

The dimensionally continued constraint equations 
EH written for the two-point-mass system read 

ViR = ( 7 ik 7 n ^ - ^-( 7 a ^) 2 ) 

+ X]( m « + 7a PaiPaj)^S a , (2.5a) 

a 

' - X<V (2.5b) 

a 

Here R denotes the space curvature of the hypersurface 
t = const, 7 b := 7 *| g (x a ) is the finite part of the inverse 
metric 7 U ('y^'Yjk = 81 ) evaluated at the particle po¬ 
sition (which can be perturbatively and unambiguously 
defined, see Appendix IA 41 of the current paper), and Dj 
is the d-dimensional covariant derivative (acting on a ten¬ 
sor density of weight one). Let us note that the source 
terms of both constraint equations are proportional to 
Dirac 6-functions. 

We employ the following ADM transverse-traceless 
(TT) coordinate conditions, resulting in irreducible 
canonical field variables, 

/ d — 2 \ 4 /f d-2 ) 

7 a = \ i + 4(d_ i) V 6ij + h7r " = 0) ( 2 - 6 ) 

where the metric function hj^ is a symmetric TT quan¬ 
tity, i.e. 


hl T = 0, djhJ? = 0, (2.7) 

and the field momentum 7 r y is split into its longitudinal 
and TT parts, respectively 

7 r7 = 7fb -|- 7 rfjtp. (2-8) 


After solving [with the usage of the coordinate con¬ 
ditions EH) and by a perturbative expansion] the con¬ 
straints EH with respect to the longitudinal variables 
4> and 7 r y , we plug these solutions (expressed in terms 
of x a , p a and hJ : T , n^ T ) into the right-hand side of Eq. 
EH- This way we get the reduced ADM Hamiltonian of 
the total matter-plus-field system, which can be written 
in the form 


H re d , p a , , ^ , TT TT ] 

= - Jd d xA d cj)[x a ,p a , ( 2 . 11 ) 


This Hamiltonian describes the evolution of the matter 
and independent gravitational field variables. The equa¬ 
tions of motion of the bodies read 


dH red 

dPa 


Pa = -- 


<7ffred 

dx„ 


( 2 . 12 ) 


and the field equations for the independent degrees of 
freedom have the form 


—/7 TT = X TTkl SH ied 
dt ij ij 5n «r : 


d _ij _ x TTij 8 H Ied 
g-^TT ~ °kl ’ 


(2.13) 

where the d-dimensional TT-projection operator is de¬ 
fined by 


S™ 3 := ^(8 ik 5ji + SuSjk) - 

2 T Sjidik T Sudjk T 8jkd%i)A d 

+ ^— -(8ijdjd + 8kidij)A d 1 

+ ^d ijkiA d 2 . (2-14) 


In the current paper we are interested only in conser¬ 
vative dynamics of the matter-plus-field system and we 
want to describe this dynamics in terms of only matter 
variables x a and p a . An autonomous (thus conservative) 
matter Hamiltonian can be obtained through the tran¬ 
sition to a Routhian description. By means of the first 
field equation (12.131) one expresses the TT part of 
the field momentum as a function of matter variables x a , 
p a and field variables h] J T , /i7 T . Then the Routhian is 
defined as 


The longitudinal part of the field momentum can be ex¬ 
pressed in terms of a vectorial function V 1 , 

n ij = diV j + djV* - jS ij d k V k , (2.9) 

and the TT part satisfies conditions analogous to EH, 

7r« T = 0 , djTr^r = 0 . ( 2 . 10 ) 


i?[x a , Pa, hJT, := H red [x a , Pa, , tt^ t ] 

— J d d xiT^hJ?. (2.15) 

Finally the matter Hamiltonian reads 

H(x a ,Pa) := i?[Xa,Pa,^ T (x a ,Pa),/lJ T (Xa,Pa)], 

(2-16) 
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where the field variables hj ^, hj 3 - are replaced by solu¬ 
tions of field equations, i.e. they are expressed in terms 
of the matter variables and eventually their time deriva¬ 
tives, which are in turn eliminated through lower-order 
equations of motion (this procedure is equivalent to per¬ 
forming a higher-order contact transformation [70117l| ). 

The dependence in the right-hand side of Eq. (12.1611 
of the field functions on the same position and momen¬ 
tum variables as those being located outside the field ex¬ 
pressions vitiates dissipation. For dissipation (radiation 
damping) to be obtained, another set of (primed, say ) 
variables has to be used in the field expressions, cf. jl4| . 
Because of our restriciton to the conservative dynam¬ 
ics, we can iteratively solve the field equations by means 
of the time-symmetric half-retarded plus half-advanced 
Green function (see Sec. IIVI below b which from the very 
beginning excludes dissipation. 

III. THE 4PN-ACCURATE REDUCED 
HAMILTONIAN 

To compute the 4PN-accurate reduced Hamiltonian 
given by Eq. (12.111) we have to perturbatively solve the 


constraint equations (12.5[) through 4PN order. To do this 
we expand these equations into the PN series, i.e., into 
series in powers in the inverse velocity of light, e := 1/c. 
We take into account that 


m a ~ 0(e 2 ), 

Pa " 

- 0(e 3 ), 


0~O(e 2 ), 

fi TT - 

n i 3 

- 0(e 4 ), 


n ij - C(e 3 ), 


- 0(e 5 ). 

(3.1) 


To compute the 4PN-accurate Hamiltonian we need to 
expand the Hamiltonian constraint equation (12.5a[) up to 
the order e 12 . It is convenient to put this equation in the 
form solved for the (flat) Laplacian of the metric func¬ 
tion (j >. After long calculation we obtain (from here the 
numbers written in subscripts within parentheses denote 
the formal orders in e) 

6 

A^ = ^ $ (2n) + 0 (e 14 ), ( 3 . 2 ) 

n—1 

where 4>( 2 ), • • ■, < f > ( 12 ) are given by 


$( 2 ) = ^E TOa<5a ’ 

a 

$ (4) = - E 
*(6) = E 


Pa - d - 2 


6a- 


4(d-l) 


*(8) = E - 


2 m a 

( p 2 ) 2 , p 2 

8 m 3 2 (d — 1 )m, 

( p 2 ) 3 ( p 2 ) 2 


> A, 


d.q>, 


S a - ( K 11 ) 2 + 

, (d + 2)p 2 a 


(3.3a) 

(3.3b) 

(3.3c) 


16to 3 4(d—l)m 3 16(d—1) 2 to q 




E 


PaiPaj 


5a~ 


d -2 


2 m a 4 (d — l) 2 


((6 - d)#,ij + 3<j) } i4>j) \hj? + -(^- T fc ) 2 - ^ h llk h Ik,j + - 2^7Ttt' 


(3.3d) 


$ 


( 10 ) 


= E 


(5(PI) 4 


+ 


3(P 2 ) 


2 ^3 


\128rn£ 16(d — l)m 3 32(d — 1) 2 to 3 


j. , (d + 6)(p 2 a ) 2 (d + 2)dp 2 3 \ (10-3d)(4-d) 2 2 

0+ 32(<i-l)^ + 96(d - !).„,/ j 18 "- 16(j-l)» * > 


+ {E(-£ 

^ a ' 

5 — d 


1 


4 m 3 (d — 1 )m, 

d -3 


A it i 3(d 2) (d 2)(6 d) 2 n ~ ik ~ik\ ,tt 

( P)PaiPaj5a+ 4 (^ _ ^3 + S(d — l) 3 ^ ^ ^ Jo 

-d 


4(d-l) 

6 — d 
2(d — 1) 


(A 


TT\2 




fprp rp^ 4 dj rp^ rprp 

'■k djj hij^k d - ^ ^ikh^j dik : j 


1 , 




TT A J.TT 

ij 


4 -d 
d- 1 


<jm' u — ( 7 rl(? T ) 2 , 


(3.3e) 


„ = V (-Ml 5 (Pq) 4 , _ 3(d + 10)(p 2 ) 3 2 _ (d + 6)(d + 2)(pir _ (d + 2)d(3d - 2)p 2 A 

^ 12 ^ \ 256m^ 32(d—l)m£ 128(d—l) 2 m 3 192(d — 1) 3 to 3 1536(d—1) 4 to 0 J a 








































d + 6 

8 {d — l) 2 ro, 


6 


(d - 3)(10 - 3d) (4 - d) 3 t ■ 2 

+ -4^1^-^ + {£( 

</> 2 Q (6 - | h. 


2\2 


+ 


3P a 


d 2 - 4 


16 m® 4(d — 1 )?ti 9 

TT / (5 — d)(10 — d) 


^ ) PaiPajda T 


d- 2 


d- 1 


~ ifc ~jk 
>7T 7T J 


32(d — l) 4 

f _ PaiPaj 


1 A f i 


7 — d 


d- 3 


2 TO a 2 (d — l ) 2 


16(d — l) 2 ar 4(d-l) 2 

30,* + ^(10 - 


(m 2 ) (hj 3 t ) 2 


TT 


10 -d ,, . 

+ 4(d-l) 2# ’ fc ^ 


-(8 — d)hj^ k — (4 — d)/i,V ) + 


TT \ (6 — d)(10 — d) 


16(d — l) 2 


3 - • TT \2 1 i TT 7 TT 


±2 I “/ill \'2 ' ill I I | 7 TT a 1.TT 


, 7i TT / ^ lTT r.TT 3 , TT .TT , . TT iTT 3 TT XT \ lTT.TTa i TT , lTT. TT / . TT ,TT \ 

U ifeAa 2 hik ’ lh i k ’ l+hik ’ lhkl ’i 4 h ^ h klj) h ik h kj d /iy +h ij h kl (h ikJl h ijM ) 


(10 3d) (4 d) _ 4 d 


) ■ 


8 (d — l) 2 ^ " " TT '‘ TT 2(d-l)^ TT; 

We also need to expand the momentum constraint equation (12.5bl) up to the order e 9 . This expansion reads 


(3.3f) 


= n ( 3) + n* (5) + ir m + nj g) + e>(e n ), 


(7) 


T9) 


(3.4) 


where I1W, ..., II) g ^ are equal to 


(9) 


^-(3) O 


2 

nI(5) = d- 1 ( 2 


n ( 7) 16(d - l) 2 


n < t > ^ ' Paidg <t>,j 
a 

d + 2 


n ( 9) 96(d — l) 3 


^ ^ ^ Pai^a 
a 

(d + 2)d o r 


d -2 ^ . 


M~ 2) 2 

16(d — l) 3 


+ jfl#* | W + * - §*&,J ) + i^rrp^^TT + ^J <T- 


L TT 


TT 


1 


h 2 d> in y — 


TT 


+ ^Pajda + & k - h^\ 

a ' ' 

( phij ^ 1 Paj da 


d — 1 


0 '‘TT’ 


(3.5a) 

(3.5b) 

(3.5c) 


d- 


d- 2 


TT 7 TT \ _jk 


(3.5d) 


We use the ADM canonical approach in an asymptotically flat space-time and we employ asymptotically 
Minkowskian coordinates. Therefore we have to assume that the functions which enter the formalism have the 
following asymptotic behavior for r —> oo (see [69[ for discussion of asymptotics in the d = 3 case): 


1 7 

r d —2 > n ij 


TT 


1 


1 


1 


r d-l ’’ 


TT r d -1 ' 


(3.6) 


Making use of the above asymptotics and the expansion m (EH) of the Laplacian A d<t>, after dropping many 
total divergences which decay fast enough at spatial infinity (so they do not contribute to the Hamiltonian), the 
4PN-accurate reduced Hamiltonian can be written as 

#<4PN [ x a, Pa, h^ T , 7T^ T ] = f ™ad a + ^ ft (2n) ( X ) X a> Pa, hj^, 7T^ T ) V (3.7) 

'' V n on —0 / 


where the Hamiltonian densities h 1 ^, ..., read 


7 red _ \ ' Pa X 
ft (4) -2^0 da 


d- 2 


2 m Q 4(d — 1) 


> Arf 


(3.8a) 
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2 (d - 1 )to 0 


*« + (* <j ) 2 , 


(3.8b) 


E) = E 


[(Pi) 3 


(P 2 ) 2 


(d+ 2 )p 2 


\16mj) 4(d — l)m 3 16(cZ — 1) 2 7 


0 2 ) S a + 


4- d 
2(d- 1) 


</>('■ 


-Jj\2 


+ -E 


PaiPaj (d-2)(d-3) 


L <5a 


0,7 ) h 


lTT 


J,red 

“( 10 ) 


-E - 


4(d — l ) 2 
3(p 2 ) 3 (d + 6)(p2)2 


iv,o i 


a 2m “ 

5(pg) 4 

128m£ 16(d—l)m^' 1 ' 32(d— l) 2 m| v " 96(d— l) 3 m a 


(3.8c) 


(d + 2)dp a 3 


2 , 3 \ , (10 — 3d)(4 — d) 2 , -. 2 

0 3 «y a + 10fJ \,o. <t> (* ) 


16(d — l ) 2 


+ 


Ef- 

< J V Zlrr 
a 

d — 3 


1 


\4m| (d — l)m Q 


PaiPaj^a 


d -2 


(d — 2)(d — 3) 
4(d — l ) 3 ‘ 

6 — d 


; o.j + 2n ik a jk \hl ; T 


4 — d 


^( 12 ) = E 


t (A^) (Z ^) 2 + ) 2 + ^07^ + ( ^ t)2) (3. 8d) 

5(p 2 ) 4 3(d + 10 )(p 2 ) 3 2 , (d + 6 )(d + 2 )(p 2 ) 2 3 (d+2)d(3d-2)p2 > 


4(d — 

fMl! 

\256m® ' 32(d-l)m7 r 1 128(d-l) 2 mf 


+ 


192(d — l) 3 m3 


(d-3)(10-3d)(4-d) 3 2 f v- 7 3(p 2 ) 2 3p 2 

- 48W _ 1)J - *■ («■ > +(U-fcr- 4(J _ 1K ^ 


1536(d — 1) 4 to q 
d + 6 


L 0 4 4 


4(d — l)rn 3 8 (d — l) 2 m a 


+ 


(d + 2)(d — 2)(d — 3) 
32(d — l ) 4 


PaiPaj x , 3(4 + 4d — d 2 ) 
O a - 


lY ’ 3 d- 1 

,2 ' 1 (6 — d )(10 — d) 


i j, l.TT.TT , (6 —d)(10 — d) 2/.TT-12 
VQm J h ik h jk d 64 (d_ X )2 ^ 


2m a " a ^ 16(d^l) 2 1 16(d—l) 2 

1 ; TT / i,TT /lTT , lTT \ i 1 iTT tTT \ , a -ifeiTT _jk , 
2 ^ij y^ikjyo-jkj + o.ji f,) +) + 4,n h^ 7r TT + 


TT - jk -*- (1 ° M)(4 E 2 7r«7r¥r + f "E (4r) 2 - 


8 (d — l ) 2 


2 (d — 1 ) 


(3.8e) 


In the next step we perform the PN expansion of the 
field functions 0 and 7 r u . From Eqs. (13.81) it follows that 
to get the 4PN-accurate Hamiltonian we need to expand 
the function 0 up to 0 (e 10 ) and the function 7 p J up to 
0 (e 9 ), 

0 = 0 ( 2 ) + 0 ( 4 ) + 0(6) + 0(8) + 0 ( 10 ) + ^(e 12 )) (3.9a) 
* ij = 4 } + 4 ) + 4 ) + 4) + Oie 11 ). (3.9b) 

Equations (13.31) imply that the functions 0( 2 ) and 0 ( 4 ) 
depend only on matter variables x a , p a , the function 
0 ( 6 ) depends on matter variables and on hj?, whereas 
the functions 0 (§) and 0 (io) depend both on x OJ p a and 
on hj^, 7 r^ T . The leading-order and next-to-leading- 
order functions 7 r^ and 7 r^ depend on matter variables 
only; the functions 7 r*^ and 7 r^ depend on both matter 
and TT field variables. The overbar in the functions 0(6), 
0 ( 8 ), 0 (io), 71(7), and 7 r^ means that they depend on the 
non-PN-expanded TT variables hj ^, 7 r^ T . 


To obtain the equations fulfilled by the functions 0( 2 ) 
up to 0( 1O ) we substitute Eqs. (13.91) into Eqs. (13.31) and 
reexpand them with respect to e. This way we first obtain 
the Poisson equations for the functions 0( 2 ) and 0 ( 4 ), 
which read 

Ad0( 2 ) = ^ ~2m. a S a , (3.10a) 

a 

(310b) 

where in the right-hand side of (13.1 Obi) we have used 
(I3.10al) . We have found it useful to split the function 
0 ( 4 ) into two pieces, 

0(4) = _ 2 ' S ' (4)1 + 4 (d- l)^ 2 ’ 


(3.11) 

















































where 5 ( 4 )! and 5 ( 4) 2 fulfill the following Poisson equa¬ 
tions 

D 2 

A d 5( 4 )i = ^ —(5 q , A d 5 (4)2 = 0 ( 2)E^. (3.12) 

a a 

In the case of the function <f>^ it is also convenient to 
split it into two pieces, 

0(6) [ x ; x a , Pa) hJj T ] = <^(6)l( x ; x a, Pa) 


+ 0(6)2 [ x ; x a,^ T ], (3.13) 

where the first piece 0 ( 6 ) i depends only on matter vari¬ 
ables x a , p a , whereas the second piece 0( 6 ) 2 depends on 
x a and (functionally) on hj^. The Poisson equation for 
the functions 0( 6 ) 4 and 0( 6 ) 2 read [in their source terms 
the Laplacians Ad4>(2) and Awere eliminated by 
means of Eqs. (13.101) ] 


a 4 V- f (Pa ) 2 , (^ + 2 )p 2 ^ (d — 2)m a ( a , d-2 /jl2 p A 

A ^(6)i - Z. + 8 (d-l)m/ (2) “ 8(rf - 1) 1 5(4)1 + 2(d^T) ^ 2 ) " ^ (4)2) J 


- (^ 3 )) 2 , 


Ad</>( 6)2 - j-0(2),i.Ay 


TT 


(3.14a) 

(3.14b) 


The Poisson equations fulfilled by the functions 0(g) and <^>(i 0 ) we present in the form where the lower-order 
Laplacians A<z0( 2 ), Ad</>( 4 ), and Ad0(6) are not replaced by their source terms. The equations read 


P 2 ( d + 2 


4-d 


A^(8) " - E ^ ^ + 4 (i Pa i)m3 </> ( 2) + 2(d-l)m a (^T)^ 2 ) “ ^ (4) J } " 2^I)^ ( ^) )2 


- 2 ^(3) f (5) 


d-2 


d-2 


l TT 


E P< 2ma 5a ~ 4(d- l) 2 ( 3 ^(2). i<?!) ( 2 )d + ( 6 “ d)<t>(2)<f>(2),ij) + ^77Y<A(4),ij| hi 

+ ^(^Sifc ) 2 - + hJ^Adh]? - 27^)71'tt - ^y( 0 ( 2 )A d 0 ( 6) + 0 ( 4 ) A d 0 (4 ) + 0( 6 )A d 0( 2 )), (3.15a) 

a i V-(5(P 2 ) 4 , 3(p 2 ) 3 , , (p 2 ) 2 M + 6 , 2 , 

A ^(io) - E j 128^1 + 16(d-l)m^ (2) + 4(d-l)m 3 l,8(d^T) *< 2 > “ ^ (4) 

Pa / (d+2)d ^3 d + 2^^ , I N \ 1 * 

+ 2(d-l)m a U8(rf-1) 20(2) " 4(d^T)^ 4) + ^ (6) J f a 


( l6(d-l) ^ 2 ) (f (3)) 2 + ^w(*(3)) 2 + ^(3)^(5) ) ~ - ^(3)\7) 


+ < 


+ 


m a \ 4m 2 d — 1 

a \ a 

6 — d / 1 


d- 2 / 3 


E Pa ! P<13 ( £72 + 73T^( 2 )) + !) _ 0(2),i0(4)b - 0(2),j0(4),i 


( 2(d — 1 ) ElEbb 0 ( 2 ),b' 0 ( 4 ) 0 ( 2 ) 0 ( 4 ),b) + 0 ( 6 ),b j 2 ^(3)T(3) [ ^ 


'TT 

bj 


4(d — 1) V2(d — 1) 

4 (d — 1 ) ( Ad< ^( 2 ))^*J )' ~ ^ _ ^(2),ijh-ik hjk ~ 2 ^ _ < t > (?),khij h i; j k + —j——^(t>(2) t khij h ik j 

d _ ^ 0 (2)7T(3) + ^(sjVtt - (^tt ) 2 


6 — d / 1 , tt z TT 3 xx 2 , TT A iTt! / 4 — C? 

4” ~ 4 t'byfc) hjj Add, 


2(d — lj w \2 lJ,fc ® fc,J 4 V o “ *j 

d-2 

4(d^T) 


(0(2) Ad0(g) + 0( 4 ) Ad<t>(fi) + 0(6) A<20( 4 ) + 0(8) A <20( 2 ))- 


(3.15b) 
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The equations satisfied by the longitudinal field momenta ^Ay ^( 5 )’ 'rfjry anc( ^(g) we obtain by substituting 
expansions (HTTTli into Eqs. (131) (1331) and reexpanding them with respect to e. The result is 


'(3) 


j — q ^ , Pa iAa.' 


# (5)J = -73T^(^( 2 ) # (3))» 


=3-3 

7 r 


(7 ),j - d i\ ~ + ^(2)^(5) 


■ 4>(2 p 


6 — d 


tt; 


21__ ~jk iTT , V k ITT 

8(d_l) 2 ^( 2) (3) ( 3 A fe + V (3) n jk,i 
6 — d 


=.n _ 

ft(4 — 


1 


(3.16a) 

(3.16b) 

(3.16c) 


(9)b = d i[ - rfTl (^(6)^(3) + ^(4)^(5) + ^(2)^(7) + ^(4)*Tr) ~ 4 (rf _ 1)2 ^(2)0(4)^) ~ ^(2)*(3) 


6 — d 


1 


g( d _“[)2 <A(2) (*(5) + 4 t) + d n^(2)^(3) /l ^ T + ^(5)^ - ^ T 4 t 

^Tl (^(2)^(3) + ( d - l )^{t)) h ^k + ^ l S4- 


(3.16d) 


r 


By making use of Eqs. (13.1 HI) (13.161) and performing where the Newtonian and the 1PN ft( 6 ) densities 
very many integrations by parts in space, it is possible depend only on matter variables x Q , p a only, the 2PN 

to rewrite the 4PN-accurate reduced Hamiltonian (13.71) in densitity ft( 8 ) depends on matter variables and on h] J T , 

the form in which its density depends on momenta p a and whereas the 3PN ft(i 0 ) and 4PN ft( 12 ) densities depend 

on the following functions: </>( 2 ), 5(4)1 > 5(4)2, V($)i ^sp on matter variables and on field variables /i^ T , n^ T . The 

^(6), hJ : T , n l ^, T [for convenience we also use and 7 r^, dependence of h(\ 2 ) on hJ 0 T is both pointwise and func¬ 

tional and this is why we have used square brackets for 
arguments of ft( 12 ). The explicit forms of the Newto 
/i( 4 ) and the lPN-level ft( 6 ) densities are as follows: 

( Pa {d~ 2 )m a 


which can be expressed by PA, and PA,, respectively, by 

means of Eq. Cm We will display now this form. ‘ The arguments of h (12) . The explicit forms of the Newtonian 

4PN-accurate Hamiltonian is the sum of pieces related 
with different PN orders, 


H red lx n h TT t r b 1 
rz <4PN L x ai Pa; a ij ; 7r TT J 

= [ d d x /i< 4 PN [x; x a , p a , hjj 1 ', 7t^ t ] , (3.17) 


1(4) 


(x;x a ,p a ) = 




ft(6) ( x ; x a i Pa) 

where 


h 

1 

^0 

h red py -x- n h TT TT*- 7 1 — 
^<4PN L x ’ x «’ Pa? 5 *• ttJ — 

7n a d a + ft( 4 ) (x; x Q , p a ) 

8(d — 

+ h(Q) (x; x a , p Q ) 

+ ft(8) ( x 5 x a,Pa,ftp T ) 

+ ^( 3 ) Pai 


\2 m a 4 (d— 

(p 2 ) 2 (d + 2)p 2 


J^( 2 ))*a» 

(3.19a) 


S(4(1 - 


8m 3 8(d — l)m Q 

(5'(4)2 - 


h 2 ) 


d-2 


+ ft( 10 ) ( x 5 X a ,Pa,^ T ,7r^ T ) 
+ ft( 12 ) [ x j x a, Pa 5 hij 1 ^Tt] > 


(3.18) 


(3.19b) 


For displaying the 2PN-level density ft( 8 ) we define two 
auxiliary functions which depend on matter variables 
only, 


k {8) ( x ; x a , p a ) := Y 


(p 2 ) 3 , (d + 2)(p 2 ) 2 


16m 3 16(d—l)m 3 


a , (d + 2)p 2 ( x 2 , 0 d-2 0 

^ 2) + 16(d-l)m a (* (2) + 5(4)1 " 2(d^l) 5(4)2 


(d-2) 2 


d-2 


3(d — 2) 


32(d - i) 2 +35(4)1 - i^ry 5(4)2 J 0(2) \ 5a ~ ^l ) )2 


(3.20a) 






























^(4)ij ( x 5 x a? Pa) •— ^ ^ 


PaiPaj c d-2 

a ^ 5a ~W^) (t){2) ’ A2) ' j ' 

By means of these functions the density /q 8 ) can be written as 

h { 8 ) (-K;x a ,p a ,hJ^) = >c ( 8 ) (x;x a ,p a ) + is , (4)ii (x; x Q , p a )/iJ T + ~(hj? k ) 2 . 

To display the 3PN density d( 10 j we define another two auxiliary functions 

( n s (d±4Kp|)^, 

(1°)( ’ a,Pa -*' 128m£ 32(d — l)m® ^ (2) 

(P 2 ) 2 


10 

(3.20b) 

(3.21) 


f(d + 6)d rk2 , fr] , d(d-2) c \ 

V 2 (d^l) 0(2) ( } (4)1 " (4)2 J 


32(d- l)m| 

(d + 2)pl f d(3d — 4) 2 , J£I (3d — 4)(d —2) \ , 

" 16(d-l)®m a V 12(d=l )' h2) (4)1 -4(d — 1) 5(4)2 J 0(2 > 

, (d—2)m a ( (d — 2) 3 ^ 4 | (3d — 4)(d — 2) , 2 (d — 2) 3 c 2 

32(d— l) 2 \4(d — 1) 2 ^ (2) + 2(d — 1) 5(4)1</>(2) " (d^lF S(4)20(2) 

2 (3d-4)(d-2)„ „ (d-2) 3 o2 \ 

+ 2(d^l) 6(4)16(4)2 + 2(d-l) 2 6(4)2 J 

y (3)P« / (3d — 2) (3d — 4) 2 , J „ (3d — 4)(d — 2) M r 

+ 2(d = l) l-8(d^l)-^ (2) + d5(4)1 -4(d — 1) 5(4)2 ) r 


3d-4 


((d- 2)5 (4)2 ,i - (3d — 2)^> (2) 0 (2) ,*) - —Y 5 (4)i,iJ V ( 3)T 


4(d — 

%)«(^.p.) : = E r i + - iv k- v m„ - + 2^3)..^),,' 


(3)’ 


(3.22a) 


4(d — 2) ■ , , (d + l)(d — 2) 


T „ , v- ■ -/r O (2d — 3)(d — 2) 2 M a 

V (3),j V (3),k + 4 (d _ 1)2 < ?(2),i A (4)l,j-g( d _ ^3-9>(2).i*(4)2,i 


(d — 2) (3d — 4) 
+ 8 (d-l) 3 

The 3PN density /i( 10 ) then reads 


0(2) 0(2),i 0(2),j- 


(3.22b) 


^(10)( x ; x a,Pa,dy T ,7r^ T ) = ^( 10 )(x;x a ,p a ) + ^ 0(2),»V(3) (0(2)703)) TT + B (6)ij( x ; x a,Pa)^ T 


1 


2 (d — 1) 


^( 2 )^ T A d ^ T - 2 ^_ ^ (2)^(3)TfTT + (^Tt) 2 ’ 


(3.23) 


where, to diminish number of terms, we have introduced the TT projection of the product 0(2)703), which, by virtue 
of Eqs. (1C9I) and (I3.16b[) . can be written as 


(0(2)703)) = 0 ( 2 ) *( 3 ) + (d - 1 )T 


*7 

(5)‘ 


(3.24) 


This TT projection should be treated as a function of matter variables only. 

To display the very large formula for the 4PN-order density h( 12 ) we introduce three auxiliary quantities 
x^i 2 ), and >r 3 12 ). We first define the function x^ 12 ) which depends only on matter variables x a , p a , 


012 ) 


(x;x a ,p a ) := ^ 


[ 7(P 2 ) 5 5(d + 6)(pg) 4 

1 256m® 256(d - l)m 7 a ^ (2) 
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(Pa) 


2 \3 


64 (d — 1 )m® 


d+1 A 


d — 1 \ 2 


-(d + 10)</>(2) — (d — 2)5(4)2 j + (d + 6)5(4)i 


2'i2 r 


+ 


(d + 6)(pg) 
32(d — 1 )m^ 

(d + 2)p 2 
32(d — l) 2 m a 


(d + 2) (3d — 2) 3 


(3d - 2)(d- 2) 


24(d — l) 2 0(2) + 2(d^I) V (d + 1)5(4)1 4(d — 1) 5(4)2 ) 0(2) + 0(6)1 


d(3d — 2) (2d — 3) 4 


48(d — 


M-3),4 , d ( toJ (d — 2) (2d — 3) c \ , 2 

l) 2 - 0(2) + 4(d — 1) [} M ~ 2 ’ S W 1 ^^1- 5 ( 4 ) 2 J ^(2) 


(d — 2)(3d — 2) 


(d — 2) 2 (2d — 3) c2 


+ 2 V + 1 ^( 4 ) 1 2 (d — 1)- S '( 4 ) 1 ' S '( 4 )2 +--*“*(4)2 J + (3d - 2)0 (2 )^( 6 )i 


(d — 2)m a I” (d-2) 4 5 (d — 2) 2 

32(d - l) 2 [ 16(d - l) 3 0(2 > 4(d - l) 2 

d- 2 /I 
2(d — 


, 5 (d-2) 2 ^ nI 5(d — 2) 2 

0(2) “ 4(d — l) 2 ^( 2d_ 3)6(4)1 “ 4(d — 1) 6(4)2 J 0(2) 

j f 1,0, 0^2 (d-2)(2d-3) 0 0 , 5(d-2) 3 o2 ^ 

1) l 2 ~ 2 ( A (4)1 5(4 )i5(4)2 + 4^2 *(4)2 ) 0(2) 


+ 


/5(d- 2) 


V 2(d — 1) ( 5(4)2 ~~ 0(2) ) ~ (3d - 2)5(4)!J 0(6) 


x - — _ 2 ,• ^ 2 

u a 


4(d — 1) 


0 ( 6)1 (^ 


(3 V 


id-2 
8(d — l) 2 


-(2d-3)4) + - (3d -2)5(4)!- 


(d — 2)(2d — 3) 
d- 1 


*5(4)2 ) 0(2) 


(^3)) 2 


(3.25) 


The function >c 2 12 -j depends on matter variables x 0 , p a and on the field function /i^ T (it contains terms linear, 
quadratic, and cubic in d^ T ). It reads 


^12)( x ; x a,Pa,^ T ) = i~E 


PaiPaj f 3(p If . (d + 4)p l 
8 m„, V 2m 4 (d—1 )to 2 (2) 


(d + 6)d 2 d + 4 


2(d — l) 2 


0 ( 2 ) 


d-i 5(4)1_ (d^TF D(4)2 i 0a 

0(2),i0(2),j 


d(d - 2) 


5 c 


d(d - 2) d-2 

4(d- l) 20(2) ’ ii0(6)1 “ 8(d — l) 3 


2d 3 fd + 2 2 o \ 3d 2 

2(d^l) V^“ 0(2) “ (d " 2)5(4)2 J + ~~2~ ‘ 5(4)1 


3d 5, 


(4)M ' 


(5d — 8)(d — 2) 


2(d — 


pj-^ 5 (4)2,i) 0(2)0(2),j | - 16(d- 1)2 + 3 ) 5(4)1,*5(4)1 J 


5(4)1,i5(4)2j + 


(3d-5)(d-2) 2 
4(d — l) 2 


(2d - l)(d- 2) 
d- 1 


- ^(3),^(i), fe - 2^3),^(3)P + 2 ( 1 - - ) ^ „• K 


5(4)2, *5(4)2, j 


2(d — 2) 
d- 1 


^(3)^(3), j 


d I V (3)j v (3),k 


0(2) + 2F(3) (V(3)j0(2),fe “ ^(3) > fc0(2),i 


+ 4(d — 1) 


^(3)^(6)J - ^3)^(5),* ^ ^3),^(5),i ^ ^(3),^(E + ( 1 " ^ 1 ^ ^ M 


0 


(3),k v (5),j f v (3) J '(5) 


,fc)] 

+ (E ^ - 4 i6( 8 /_ 1^2 0(2)40(2),! - I6(d _ 2 1)2 (^C 2 )^(2).« + 5 W 2 .«) + i ^ L T y5(4)! iij )^ T ^ T 


E 


3Pa 


-da 


20-d 2 

16(d — l)m a a 64(<i— 

d 4 \ mnn nrnr 1 


4(d — 1) 

^ / a \2 N \/. TT)2 1 f (d + 2)(10 - d) 2 

^(0(2), fe ) JK- ) 8(d _l) ^(2)+ 35 (4)l 


4 (d_ i) 5 (4)2 j^ T Arf/ig T - -(2 hj^ (dJ fc T , + hjj^) + h^h^hJT. 


(3.26) 


The function >£? 12 , is proportional to the second part 0(6)2 of the function 0(6). It is a function of matter variables 

nd, through the function 

(d + 6)(p 2 ) 2 , (d + 2)(3d — 2)p 2 


and it depends on /id 1 both pointwisely and, through the function 0(6)2, functionally. It equals 


*fi2)[ x ; x a,Pa,dy T ] : = S E 


32(d— l)m3 32(d— l) 2 m a 


~0(2) 
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, (d - 2 )m a {5(d - 2 ) 2 (o x2 ^ fnJ ONO 

+ 32(rf- 1)2 U(d-l) {S{4)2 ~ 0 < 2 >' (M “ 2)5(4)1 

The 4PN-order Hamiltonian density d( 12 ) finally equals 


Sa ~ 4 (d- 1) ^(3)) 2 + 4 (rf_ 1 )2^(2),ij /l 5 T |^ > (6)2- (3-27) 


^(12) [ x i x ai Pai hij i ^Tt] — ^( 12 ) ( X i X ai Pa) + ^ 1 2 ) ( x > x ai Pa> ) ~b ^(12) [ X ! X a > Pa, h^j ] 

, d(lld — 18) , 2 0 (5d — 8)(d — 2) \ ,, _ij >tt 3d— 5 _ %J . 

+ V 8(d- l) 3 0 < 2 > + d^l 5(4)1 4(d — l) 3 5 ( 4 ) 2 J 7r (3)W(2)^(3)) - (dTi)3^(2) (^(2)^(3)) 

d — 3 


ij \TT 


+ 2 ( 2 ^(3) /l 5 T + ^)( 2 ^ “ *©) (^TT - |( < / > (2)^(3)) TT 


d- 1 


0(2) Kt) 2 


/ (9d — 14) (4 — d) 
V 4(d — l) 2 


0(2) 0(2) ,j V (3) ~ JZTJ^ 4 ) 1 >J V (3) + 


(3d — 4)(d — 2) 


2(3d — 5) 


4(d_l)2 3(4)2,jV( 3) + d l 0(2), j ^(5) 

(3.28) 


IV. FIELD EQUATIONS 

Dynamical degrees of freedom of gravitational field de¬ 
scribed by the functions hj^ and are solutions of 
the field equations (12.131) . From the 4PN-accurate re¬ 
duced Hamiltonian given in Eq. (13.171) one can derive 
4PN-accurate approximate field equations. In the rest of 
the paper we will only need to use field equations which 
follow from the 3PN-accurate part of the Hamiltonian 
(13.171) . It reads 

H red Tx r> /) TT 77 -b 1 

•“^PN L Xa ’ Pai'Hj ) ^ttJ 

= J d d a;/i<3PN(x;x a ,p 0 ,d^ T ,7r^ T ), (4.1) 

where 

^< 3 PN (^5 ^a 1 Pa i h^j , ^T r p r p) — ^ ^ ^Tla^a ^(4) (^5 ^a i Pa) 

a 

+ h(fi) (x; x a , p Q ) + d (8) (x; x a , p a , hj^) 

+ tylO) ( x ; x ajPa, hJT, 7T^ t ). (4.2) 


For this Hamiltonian the field equations (12.131) take the 
form 


ITT _ rTTfcZ <9^<3PN 

13 ~ 0i3 c>4V 


+ 0(e 7 ), 


(4.3a) 


‘TT 


,-TTj.j f 0h<3PN f <%<3PN \ 

\w 

+ ( dh TT N ) )+ C ’( e8 )> ( 4 -3b) 

V U L kLmn / ,mn ) 


or, more explicitly, 


lTT _ r 

n ij — °ij 


TTkl 


27Txt - 


2(d — 2) 
d- 1 


0 ( 2 ) 


-fcZ 

r (3) 


^(e 7 ), 


(4.4a) 


* 

t tt 


= -5. 


+ 


TTkl 


| 2^(4)fcZ ~ 2^ d ^ T + B (6)kl 


1 


2(d — 1) 

+ 0(e 8 ). 


(t>(2)^dhk7 + 0( 2)h 


, TT 


(4.4b) 


By combining these two equations one gets the equation 
fulfilled by the function hj^. It can be written in the 
form of the wave equation, 

n d+1 tif? = Sl ; T , (4.5) 

where Dd+i is d’Alembertian in (d + l)-dimensional flat 
space-time, 




:= -5 t 2 


J d+1 •- 

and where the source term is 


(4.6) 


Qk t _ r 
1 3 


TTkl 


S (m + 2.B(6)fcZ 
+ c i_i ^0(2 )AdhJ? + A d ^( 2 )/i 
+ p<9*(0( 2 )^) | +0(e 8 ). 


TT 

kl 


(4.7) 


After solving field equation (14.51) for tif^ one can obtain 
the TT field momentum 7r^ T from Eq. (I4.4al) : 

1 j t"t“ d 2 

2 


‘TT 


.jT T + £_| { ™( 0(2 ) i g>) +O (P). (48) 


In the current paper we are interested only in conserva¬ 
tive dynamics of the matter-plus-field system; therefore, 
we use a time-symmetric (half-retarded half-advanced) 
formal solution of the wave equation (H3i We then ex¬ 
pand this solution with respect to the retardation. This 
way we obtain the near-zone PN expansion of h ]’- T of the 
form 

^ = (A-/ + Afdl + A - 3 5 4 + ■ ■ ■) S^. (4.9) 
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By making the expansion (14.91) we exclude from the 4PN- 
level near-zone metric the nonlocal-in-time contribution 
coming from tail effects (what was found in Ref. [fiUp. 
Expansion of /i.T T to the order required in our computa¬ 
tions reads 

hi? = hJl 3 +hfl 3 +0{e*)-, (4.10) 

we thus omit the non-time-symmetric parts and 

. After plugging the expansion (14.101) into Eq. (14.51) 
one obtains equations fulfilled by the functions hj^ and 
hjtyij- The leading-order function is the solution 

of the Poisson equation 

= (4-ii) 

where the source term S^j is defined in Eq. (I3.20bl) . 
The next-to-leading-order function fulfills equation 

+ ^? W , (4-12) 

where the source function S(p)ij equals 

S( 6 )ij = 2 B (6)ij + j—j ^ (2 ) A d hJ^ i:j + A,(^ (2) ^,)) 

+ 2 JZJ 9t (^(2)^(3)) > ( 413 ) 

and the source function B^j is defined in Eq. (I3.22bl) . 

V. 4PN-ACCURATE ROUTHIAN 

Introducing Routhian description is an intermediate 
step on our route to the Hamiltonian which depends only 
on matter variables. The 4PN-accurate Routhian i?< 4 PN 
we obtain from the 4PN-accurate reduced Hamiltonian 
H<. 4 pn °1 Ed- HR-171) by performing the Legendre trans¬ 
formation with respect to the field momentum 7 t^ t [see 
Eq. (12.151) ]. The result is 

.R<4PN [ x a, Pa , , LJj L } = #< 4 % [x a , p a , hj ^, 7T^ T ] 

— J d d xhf^Tili T , (5.1) 

where on the right-hand side the field momentum n^ T 
is expressed in terms of /i7 T , hJ 0 T , and matter variables. 
The Legendre transformation is to leading order realized 
by Eq. (l4~8l) . We will show now that this leading-order 
formula is enough to get the 4PN-accurate Routhian. 

Let us split the field momentum in the following 
way [cf. Eq. (14.81) ]: 

4t = \ (42)43) ) TT + Sn TT> ( 5 ' 2a ) 

= 0(e 7 ), 5 tt^ t = 0(r l ~ d ) for r —> oo. (5.2b) 


We will show that the density of the 4PN-accurate 
Routhian ED does not depend on The only 

part of the density which could depend on 8tt^ t [see Eq. 
(13.231) ] consists of three terms, 

* : = (^ t ) 2 - ^ ^( 2 )^( 3 ) n TT - hJ^TT- ( 5 - 3 ) 

Making use of the splitting (I5.2al) and the relation (13.241) 
we rewrite the quantity <5c in the form 

* = -J(^7) 2 - ^3T^(2)43)^ T 

+ (57r^ T ) 2 + 2 (d — 

+ (jz l) (42)^(3)) TT (( d - - 42W(3)) • 

(5.4) 

By virtue of the representation (ED we rewrite Eq. El 
as 

6v = -\( h ¥) 2 - 

~(jZ l) ^(2)^(3) (^(2)43)) TT 
+ (<5 tt^ t ) 2 + (9 fc ED 3fe , (5.5) 

with 

ED 3fe := V(g) ^4(d — 2 )Sn^ T + 2^—-^- (0(2)^(3)) ) ■ 

(5.6) 

Because ED 3 *, decays as l/r 2d 3 for r —> oo [the quantity 
VA = 0(r 2 ~ d ) when r —>■ oo], it does not contribute to 

the Routhian, and because (^7r^ T )“ is of the order of e 14 , 
only the first three terms on the right-hand side of (15.51) 
contribute to the 4PN-accurate Routhian. 

By virtue of the above result, to get the 4PN-accurate 
Routhian it is enough to eliminate from the 4PN-accurate 
reduced Hamiltonian (13.171) the TT part 7r(j? T of the field 
momentum by means of the relation (14.81) . Making use 
of Eqs. (13.191) (13.281) and Eq. (14.81) . the 4PN-accurate 
Routhian reads 

-R<4PN [x a : Pa ■ ^ ] 

A d x t< 4 PN [x; x a , p a , /iT T , /iT T ], (5.7) 

where 

r<4PN [x; x a , p a , hJ 0 T , hj^] = ^ m a 8 a + h (4) (x; x Q , p a ) 

a 

+ h (6) (x; x a , p a ) + /i (8) (x; x a , p a , hj^) 

+ t(io) (x; x a , p a , /iT T , /iT T ) 
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+ r (i 2 ) [ x ; x a , p a , hj?, hij T ] ■ (5.8) 

The Routhian densities /i( 4 ), h( 6 ), and h ( 8 ) are identical 
with the corresponding densities of the reduced Hamilto¬ 
nian (13.171) and they are given in Eqs. (13.191) and (13.211) . 
The 3PN density t( 10 ), after some more integrations by 
parts, can be written as 

t(10)( x ;X a ,P a ,hJ?,hJ?) = *C(10) (x; X a , p a ) 

+ 20(2), iV ( J 3) (</>(2)7l(3 ) ) TT + £ (6) ij(x; x a , Pa )^ T 
+ 2(d-l) * (a)h V TAdh * F + 

-J(^S T ) 2 - ( 5 - 9 ) 

The 4PN-level Routhian density t( 12 ) can be easily ob¬ 
tained by replacing in the 4PN-level Hamiltonian density 
h( 12 ) [given in Eq. (13.281) ] the field momentum 7 r^ T by the 
two first terms from the right-hand side of Eq. (I5.2a,l) , 


We first split hj^ into the leading-order term hj^. 
and the rest Shl T : 

L J 

h]? = hj? )ij+ 5 hj?, 

= 0(e 6 ), ShJ? = 0(r 2 ~ d ) for r -> 00 . (6.2) 


The part of the 2PN density h ( S ) which depends on hj^ 
reads [see Eq. (13.211) ] 

<^1 : = )i>y T + \(hj,lk) 2 , (6-3) 


where S^ij is the source term for [see Eq. (gUJ]. 

Making use of the splitting (16.21) . Eq. (16.31) can be rewrit¬ 
ten as 



+ ^® 2 + ^%)^ T (6-4) 


r (12) [x;x a ,p a ,hJ?,K, 


W] 


We still rewrite Eq. (El) in the following way: 


= h 


( 12 ) 


x; x a , p a 


■ tt 


-h TT 
2 lJ 


d -2 
d- 1 



VI. 


(5.10) 


4PN-ACCURATE CONSERVATIVE 
MATTER HAMILTONIAN 


Sh i = \S{A)i 3 h T J )i3 + + i(A dhj^ ij )6hll r 

+ \ h jA)ij,k 6h Jj^k + \( Sh Ij^k) 2 

+ ^(5(4 (6.5) 


In this section we consider the 4PN-accurate conserva¬ 
tive matter Hamiltonian, which depends only on the par¬ 
ticle variables x a , p a . To obtain the conservative Hamil¬ 
tonian, from now on we use only the time-symmetric part 
of the field function hj ^. To get the matter Hamiltonian 

we eliminate hj^ (and hJ J T ) by replacing it by (time- 
symmetric) solution of the field equation (14.511 . We thus 
treat hj^ (and its time derivatives) as functions of mat¬ 
ter variables x a , p a . All time derivatives of x a , p a in¬ 
volved in h] 1 T (and /i? T ) are eliminated by means of the 
lower-order equations of motion. 

We thus start from plugging into the 4PN-accurate 
Routhian (15.71) the time-symmetric solution of the field 
equation for the function hj^. This way we obtain the 
4PN-accurate conservative matter Hamiltonian 

rr matter\ 

IJ <4PN P a) 

•— -^<4PN [^Xa? Pa? h%j (x, X a , Pa)? (^, X a , Pa)] ? 

( 6 . 1 ) 

where we have assumed that all time derivatives of x a and 
p a present in hjj 1 and hj? were eliminated by means of 
lower-order equations of motion (to get the 4PN-accurate 
results we only need to use Newtonian and 1PN equations 
of motion). We reorganize now density of the Hamilto¬ 
nian ED by employing the field equation ED for the 
function hJl. 

l J 


where we have employed the field equation (14.111) for the 
function hj^. Finally, making use of the identity 

(AdhJ^Shl^ = d k (hll )ijk 5hl^)-hll )ijk 5hl^ ( 6 . 6 ) 

equation (16.51) takes the form 

bhi = + -(ShJHk) 2 + 9fcEDi fc , 

(6.7) 

where the exact divergence EDit is defined as [we have 
used here the explicit form (12.141) of the TT projection 
operator] 

ED lfe := jtk 6hJ^ + Shl k T [A^S imj 

~ 2{d-l) AdlS(4)jj ’ i ~ 2(E) Ad2S(4)j,j,i )' 

( 6 . 8 ) 

Because EDi*, decays as l/r 2d ~ 3 for r —> 00 , it does not 
contribute to the Hamiltonian. Therefore the Hamilto¬ 
nian density 5h\ can be replaced by 5h\^ where 

6h[ := + \( Sh T^k) 2 , 

(Shl^f = 0(e 12 ). (6.9) 
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Because 

(*0 2 = Sh^(A d Sh^) + d k ET) 2k , 

ED 2fe := 5hJ?5hJ? k , (6.10) 

and ED 2fc = 0(r 3 ~ 2d ) for r —> oo, instead of the density 
HD one can also use 

Sh'l := \s {A)i3 hJ 4 (AdShJ*). 

( 6 . 11 ) 

The part of the 3PN Routhian density tio [see Eq. 
(15.91) ] which depends on hj? or hJ 0 T reads 

+ 2(^1) ^)^^ 

+ ~ j(^ T ) 2 - ( 6 - 12 ) 

Let us denote the sum of these terms and of the expres¬ 
sion (16.111) by Sh 2 , 

Sh 2 := ^S (4)ij hf 4 T }ij + ^(/ijj'y.fe) 2 + B (e)ijhJ^ 

+ 2 (rf - i^ ) h ^ Adh ^ ~ \ 6h l?^d 6h l?) 

+ - j(^ T ) 2 - ( 6 - 13 ) 

In the held equation (14.51) we split the source term Sj^ 
into the leading-order part and the rest SS]j T , 

a d+1 hj? = . + SSjj T , SS^ T = 0(e 6 ). (6.14) 

After substituting the splitting EH into (16.141) we obtain 
(remembering that D^+i = — d 2 + A d ) 

A d hJ? = h?? + + SS ^ T . (6.15) 


With the aid of the above relation, after dropping the 
total time derivative on the right-hand side of Eq. (16.171) , 
the density Sh 2 from (16.131) can be replaced by 

Sh' 2 := ^S {4)ij hJ 4 T )i0 + ^(hj? )ij k ) 2 + B (6)ij /i?- T 

+ 2 (d- i^( 2 ) h ^ Adh ^ 

+ 2 rf_ f ^(2),^(3)^ T - 

+ - ^5h^8Sl T . (6.18) 

The expression (16.181) contains all terms related with 
the 2PN-level and 3PN-level parts of the Routhian (15.71) 
which depend on hj^ or hj^ . Taking into account this 
expression and all other terms entering the Routhian 
ED, one can write the 4PN-accurate conservative mat¬ 
ter Hamiltonian ED in the form 

ff<4PN r (Xa, Pa) = J d d x /i^ r (x; x a , p a ), (6.19) 
where 

h< 4 pN r (x;x a , Pa ) = y m a 5 a + h (4) (x; x a , p a ) 

a 

+ *-(6) ( x ; x a. Pa) + ^““(x; X a ,p a ) 

+ /l(io) ter (x; Xa, Pa) + ^(l 2 ) ter (x; Xa, p„). (6.20) 

Here the densities /i( 4 ) (at the Newtonian order) and 
(at the 1PN order) can be found in Eqs. (13.191) and the 
densities at the 2PN, 3PN, and 4PN orders are as follows: 

^(8) tter ( x ;X a ,Pa) = >f(8)(x;x a ,p a ) + ^(4)^4- 

+ (6.21a) 


By virtue of Eqs. (16.21) and (14.111) . Eq. (16.151) finally leads 
to 

A d 8hJ^ = hj? + 5S^ j T . (6.16) 

Making use of Eqs. (16.21) and (16.161) . we can prove the 
following relation: 

-^ShJ^{A d 5hJ^) 

- - j' 5ft « T,5S 5 T - 

(6.17) 


^(io) ( x ; x . 


-a i Pa 


~ ii \ TT TT 

0( 2)7T( 3 )) 

2) 4 T A d 4 T 


( 6 - 2ib ) 


i 


• TT d 2 , _y \ TT 

2 h ij +^7lW 2 ) ,r ( 3)) 


Shf^SS^. 


’(3) 
(6.21c) 
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VII. 4PN-ACCURATE NEAR-ZONE 
CONSERVATIVE MATTER HAMILTONIAN 

We employ now the crucial near-zone (time- 
symmetric) PN expansion (14.101) of the field variable 
hjj 1 ', i.e., we use hj .) T = + ShJ 0 T with ShJ? = 

+ 0(e 8 ). We also employ [see Eq. (16.141) ] 


Here again the densities /i( 4 ) and h ( 6 ) can be found 
in Eqs. (14.191) . the 2PN density /i“™ ne (x; x a , p a ) = 
/i^ tter ( X ; x Q , p a ), where h^ tter is given in Eq. (16.211) . 
and the 3PN/4PN densities h^- zone and h^- zone fol¬ 
low from the sum of densities h™ a * ter + after plug¬ 

ging the expansions (14.101) and (17.11) into the latter. The 
3PN-level density reads 


xeTT _ oTT 
0b ij ~ *(6 )ij 


0(e 8 ), 


(7.1) 


_ near-zone 


‘( 10 ) 


(x; x a , Pa) = ^( 10 ) (x; x a , p Q ) 


where the next-to-leading-order source function S'(6)ij is 
defined in Eq. (14.131) . After plugging (14.101) and (17.11) 
into the 4PN-accurate conservative matter Hamiltonian 
of Eq. (16.191) . we obtain the 4PN-accurate near-zone con¬ 
servative matter Hamiltonian, 


rrnear-zone / 
11 <4PN 


Pa) 


J d d x h 


near-zone 

<4PN 


(x;x a , 


Pa), 


(7.2) 


+ 20(2), »Vj;3) (0(2)71(3)) + B (G)ij h J*)ij 

+ 2(d - l)^ (2)h ^ ijAdh ^ ij 

+ - l([[/74]] 0 )t 

(7.4) 


where 

h< e 4PN° ne ( x ; X 0 ,p 0 ) = ^2 m a&a + h(4)( x ; X Q , p a ) 
a 

+ 6) ( x ; X a , Pa) + /i(8) arzone (x; x a , p a ) 

+ Xa, Pa) + /l“ ne (x; Xa, Pa ). 

(7.3) 

I 


where the notation [ (Ipj ] ] 0 means elimination of time 
derivatives x a , p a in by means of Newtonian equa¬ 

tions of motion. 

The 4PN-order density we split into two parts, 

/l(lff ZOne ( x ; x a, Pa) = h\ 12 ) ( x ; Xa, Pa) + h\ v2) (x; X a ,p a ). 

(7.5) 

The first part reads 


fr( 12 )( X ; X a,Pa) = ^( 12 )(x;Xa,Pa) + *fl2)( X ! X a, Pa, + ^ 12 ) [x; X a , p a , 


0 ( 2 ) (( 0 ( 2 ) 71 ( 3 )) TT ) 


1 ((3d — 2)(3d — 4) 

d^i 1, 8 (TT) 0 < 2 > + dS{4)1 4 (d 1) 


(3d-4)(d-2) \ , _ij \tt 

S { 4)2 ) ^(3) (0(2)71(3)) 


3(d 2) 2 7 (d 2) 2 \ ^ ■ 4 k TT 

4(d- 1)2 ^( 4 ) 2 ,1 ~ A (4)1j + 4 ( rf _ 1} 2 0( 2 ) 0(2),1 ) * / (3) + 2 (« - 2)0(2),2 1/(5) - - P( 3 ), fe h(4)ij 


+ 2^) ifc ^ ife - [[^-]] 0 - ^~~jj 0(2) ( [ [^(4fa- ] ] o) 2 + 7170(2)^3) [ffl]] 



(7.6) 


The density h| 12 , depends on the following functions: 
0(2), 5(4)1 and 5(4) 2 [they determine, via Eq. (13.111) . the 
function 0(4)], VW and VW [they define, through Eq. 
( 12131) , the functions nL and 7fH>], and 0( 6 )i, 0(6)2- 
The function 0(6) i is the solution of Eq. (13.14al) : the equa¬ 
tion fulfilled by the function 0(6)2 one obtains from Eq. 
(I3.14bl) for the function 0(6)2 after replacing in the latter 
hj? by the leading-order hj^, 

Ad0(6)2 = d _ 1 0( 2 ) ,ijhj^ ij . (7.7) 


Formulas needed to express in an explicit way all these 
functions in terms of x a , p a , and x are given in Appendix 
□ The last two terms in Eq. (USD are related with the last 
two terms in Eq. (17.41) : 1 means here elimination 

of the time derivatives x a , p a by means of 1PN equations 
of motion. In the rest of this section and in Sec. IVIlXI all 
time derivatives of x a , p a are eliminated by means of 
Newtonian equations of motion; therefore from now on 
we drop the notation [[ • ]] 0 - 

The second part h 2 12 j of the 4PN density h"® 2 (' zone 
is the 4PN-level contribution coming from the density 
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(I6.18|) : it equals 


where 


^( 12 )( x ; x a , Pa) = fiwu 

+ 2(d-l) ^ 2) Adh mi + ^ A ^W>J 

+ (-B(6)*j + ^ ( 7 - 8 ) 


( 7 - 10 ) 

After plugging Eq. (17.91) into Eq. (17.81) and shifting some 
time derivatives!^ the density /i 2 |2 , can be rewritten as 


,. t 2 2 1 2 2 2 3 

The function according to Eq. (14.121) . is the sum /i( 12 ) = /i^ 2 ) + ^(12) + ^(12)’ (7.11) 

of the two terms, 


J.TT _ ^TT , a-1VTT 

ft (6)d - °(6)ij + A d ft (4)b> 


(7.9) whertfl 
I 


t.2.1 _ 4 

( 12 ) “ 2(d-l) 


(7.12a) 


, 2,2 _ 

rt (12) — 


1 


d- 2 


2( d-i)^ ) ^w« - + { 2 ^^^ + 2zrm*>*%)) ) A d lh m 

+ 7?(6)ijA d (7.12b) 


; ? 2 ’ 3 -I|V TT /) TT ? TT (A ~ 1 h TT 'A 

ft (i2) - i;b “ *(6)y( A d ft (4)ijdJ- 


(7.12c) 


VIII. COMPUTATION OF THE 
4PN-ACCURATE NEAR-ZONE CONSERVATIVE 
MATTER HAMILTONIAN 

The bulk of computations we did to derive the explicit 
form of the 4PN-accurate near-zone conservative mat¬ 
ter Hamiltonian, i.e. to perform integration in Eq. 113, 
was performed in d = 3 dimensions, where our working 
horse was Riesz-implemented Hadamard regularization 
supplemented by a Hadamard “partie finie” concept of 
a function at its singular point. The results of global 
(i.e., extending to the whole R 3 space) three-dimensional 
integrations were then corrected in two different ways: 
(i) the UV divergences were locally (i.e., within small 
balls surrounding particles’ positions) recomputed by us¬ 
ing dimensional regularization; (ii) the IR divergences 
were also “locally” (i.e., outside a large ball enclosing 
particles, or in a neighborhood of r = 00) recomputed 
by introducing an additional regularization factor ( r/s) B 


1 Shifting time derivatives means replacing AB by —AB. This 
implies adding a total time derivative to the Hamiltonian density 
and on the level of Hamiltonian is equivalent to performing a 
canonical transformation. 

2 The density , after ignoring its last term, is identical to 

the density r| pN introduced in Eq. (3.4) of Ref. l5lll (the reason 
of omission of this term in ffd is explained at the end of Sec. 
IV1IIC2II . 


(with a new IR length scale s), which modifies the behav¬ 
ior of the part of the field function which diverges 

at r —> 00 [see Eq. (IA40I) ]. The details of regularization 
procedure are explained in Appendix [A] 

Dimensional regularization introduces a natural length 
scale £q, which relates gravitational constants in d and 3 
dimensions [see Eq. (IA26I) ]. As explained below correc¬ 
tion of the UV divergences by means of dimensional regu¬ 
larization produces some poles proportional to l/(d — 3) 
together with f 0 -dependent logarithms ln(n 2 /£o). All 
these poles and logarithms were removed by adding to 
the Hamiltonian a total time derivative. The terms which 
contribute to IR divergences depend, after regularization, 
on the IR length scale s. Moreover, as we also explain be¬ 
low, different ways of regularizing IR divergences lead to 
different results, so the final result of IR regularization 
is ambiguous. We were able to express this ambiguity 
in terms of only one dimensionless constant, which we 
denote by C. 

The UV divergences is not the only problem caused by 
usage of distributional sources. General relativity uses 
standard algebraic and differential calculus of ordinary 
functions, which includes e.g. the Leibniz rule. By intro¬ 
ducing distributional sources we violate this framework: 
in some places we have to use distributional differentia¬ 
tion which does not fulfill the Leibniz rule [see Appendix 
IA 5| . The field equations (13.101) . (13.121) . and (13.141) (13.161) 
tell us that for the second (space and time) derivatives of 
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</>(„) and for the first (space and time) derivatives of ItVn 
the distributional derivative has to be applied. The same 
also holds for and 7T/^ XX , respectively, taking into 

account the Eqs. (14.111) (14.121) and (14.81) . The differen¬ 
tiation of any product of these functions goes with the 
Leibniz rule. 

The 4PN-accurate near-zone Hamiltonian is the sum 
of Hamiltonians at different PN orders, 

#< 4PN ( x ai Pa) C) = ) y TO a + H N (x a , p Q ) 

a 

+LflP N (x Q , p Q ) + i?2PN(x a , Pa) 

+7J 3 p N (Xa,Pa) +i? 4 pN’ Z ° ne (S) (x a , p a | C) , (8.1) 

where we have introduced notation which indicates 
that the 4PN near-zone Hamiltonian depends on the 
IR regularization scale s and on one dimensionless 
constant C parametrizing ambiguity in the regular¬ 
ization of IR divergences. Hamiltonians through 
H 3 pn were uniquely recomputed by us using three- 
dimensional Riesz-implemented Hadamard regulariza¬ 
tion supplemented by dimensional regularization. These 
Hamiltonians do not develop IR divergences. 


A. Computation of the 3PN-accurate Hamiltonian 


The three-dimensional 3PN-accurate Hamiltonian, 


#<3PN(x Q ,Pa) = ^ m a + 4L N (x a , Pa) 
a 

+ iLlPN (x a , Pa ) + ^2PN (x a , Pa ) 

+ H 3 p N (x a ,Pa), (8.2) 

was computed for the first time in the series of papers 0- 
HHj. We have recomputed it here and got the result iden¬ 
tical with the previously obtained. The explicit formu¬ 
las for the three-dimensional Hamiltonians ILn through 
H 3 pn written in general reference frame are given at the 
end of the current section. 

In the computation of the 3PN/4PN Hamiltonians we 
eliminate time derivatives of x a and p a [present in Eqs. 
EH and EH)] by means of lower-order equations of 
motion: we use Newtonian [in both eh and EH] or 
1PN [ in EH] equations of motion. To perform dimen¬ 
sional regularization of UV divergences we have to use 
the d-dimensional version of these equations, which fol¬ 
lows from d-dimensional Hamiltonians. We have explic¬ 
itly computed the Newtonian and IPN Hamiltonians in 
d dimensions. They read [notation “+(1 -n- 2)” used here 
means adding to each term another term obtained by the 
exchange of the bodies’ labels] 


id N (x a ,p a ) 


P? 


2mi 


n(d — 2) TO 1 TO 2 

4(d-l) 


+ (1 fA 2), 

(8.3a) 


ddlPN (x a , Pa) = - 


(PI) 


2\2 


8 m\ 4(d — 1) \2 


-(3d - 2)(pi • p 2 ) 


- d —p\ + hd - 2) 2 (ni2 • Pi)(ni2 • p 2 ) 
m 1 z 

I ^{d-2) 2 mlm 2 
+ 8(d — l) 2 r 2 2 _4 +( 

where we have introduced 

r(d/ 2 -l) 


d -2 


' 12 


k := 


47T^/ 2 


(8.3b) 


(8.4) 


Let us note that k = 1/(47t) in d = 3 dimensions. 


B. Integral of d( 12 ) 


The part h} V2j of the 4PN density given in Eq. (17.61) is 
UV divergent but it does not develop IR poles with the 
exception of the term oc 0(2) (^( 4 )ij) 2 - This term contains 
however the multiplication factor d — 3 which causes the 
IR divergence of 0(2) (^yp^) 2 to not contribute to the 
Hamiltonian (but the UV divergence of this term pro¬ 
duces some nonzero contribution). After employing the 
regularization procedures described in Appendix [A] we 
have obtained 


r. 1 

112) 


H-% 1 (x a , p a ; d) = J & d x h\ 

= X(i 2 )( x a:Pa) + X(i 2 )( x a,Pa)ln 
X(l 2 ) ( x a; Pa) 


ri2 

£0 


d- 3 


+ 0(d- 3), (8.5) 


where for convenience we have introduced a new 
dimensional-regularization (DR) scale Iq related with the 
original scale £q defined in Eq. (1A26I) by the relation 

to = v-^=e _7E / 2 . (8.6) 

ZyJ 7T 

Here 7 e denotes the Euler-Mascheroni constant. 


C. Integral of d 2 12 ) 

To compute the integral of d 2 12 ^ we have to study 
asymptotics of the function hj^ for r —> 00 . The func¬ 
tion hj^ is the sum of two parts, and 

[see Eq. EH]- By virtue of Eq. (14.131) the function 
[defined in Eq. (17.101) ] can be written more explicitly as 

c lJij = (jZT^ 2 ) h V)ij + 2( rf _ 1 dt A d 1 (^(2)*fi)) 
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1 \ 11 

+ 2A d 1 B( 6 ) ii + — —A^ (<j>(2)AdhJ^ i: j) J 

(8.7) 

Analysis of this equation shows that = 0(r 2 ~ d ) for 

r —> oo. It is also not difficult to check that the inverse 
Laplacian = 0(r 4_d ) when r —> oo (i.e., it 

diverges linearly in d = 3 dimensions). This behavior of 
C££ . and A^hJ^j for r —> oo implies that the term 
h 2 (i 2 ) is convergent at infinity, whereas the terms h^ 2 ) 
and h 2 ([ 2 ) develop IR divergences. 

I 


1. Integral of h 2 /^ 

The integral of is convergent at spatial infin¬ 

ity but it is difficult to express its integrand as an ex¬ 
plicit function of x (and of x a , p a ). To do this we 
further transform the density h^ V2 y after employing Eq. 
(ill), by making some additional integrations by parts 
in space and shifting time derivatives. We also replace 
Ad{(f > ( 2 )hJ^ i j) TT [which comes from the first term of Eq. 
(17.1 2 a.f) ] by more explicit expression using the definition 
(12.141) of the TT projection operator. We finally obtain 


d -2 


*■(12) - “( d _ X )2 (0(2)7T(3)) TT + 

d -2 


- (0(2),fe) 2 (*.(4)’ i j) 2 

(d- 2) 2 


+ JZJ <t>V ).« A d 1 (^(2 )M h mi)) - ^(2)^(3) 9 ‘ 2 ( A d 1 (^(2)^(3))) 


TT 


2(d- 2) 


( B ' 


d- 1 

where the last term equals 


( 6 ) 1 ? 


2(d- l) 


^(2) A * h m) 9 1 ( A d 1 (^(2)^(3))) + /l ?i 1 2) 1 


( 8 . 8 ) 


^(iz) 1 “ JZTJ (%)b + 2 (d- 1) ^ (2) Ad/l Mb) 

+ (-B( 6 )ij + ^ 0(2) A d^4Xj) ( A d 1 ( S (6)ii + 0(2) A d /l Wij)) 


TT 


(8.9) 


i d 7 2,2 — 'Z.'Z'L ( /^i\ 

d xh^ — X(i2)amb( Xa ’ P a ’ W 


2,2,3/ \ 

7*12 X(i2) 


. 2,2,2 / \ i 

+ X(j 2 ) ( X Q,Pa)ln 


I- 

It is crucial to single out the last term in Eq. (18.81) and to ing use of the procedures described in Appendix [A] one 

put it exactly in the form of Eq. (18.91) . Due to this it is gets 

possible to compute fully explicitly (in d = 3 dimensions) 
all inverse Laplacians involved in Eq. m • Formulas 
needed to calculate in d = 3 all inverse Laplacians and 
perform TT projections involved in Eqs. m and EH 
are given in Appendix [C] 

After performing regularization by means of the pro¬ 
cedures described in Appendix [A] we have 

*4PN 2,1 ( X a>Pa; d ) = J ^ X ^{12) 

= xll 2 )i x a, Pa) + xli 2 )( X a,Pa) In 

2,1,3/ \ 

X(ia) W-Pa) 


d -3 


+ X(i2) 4 ( x “>P«) ln ~ +0(d-3). (8.11) 

Because of the ambiguity of the results of IR reg¬ 
ularization discussed in the Appendix IA 31 the term 
X(i 2 )amb( Xq ’ Pa) C) depends on some undetermined di¬ 
mensionless constant C which parametrizes this ambi- 
guity. We have found that the coefficient X(i 2 ) °f the 
logarithm ln(ri 2 /s) can be written, after adding a total 
time derivative, in a very specific form. Thus we have 
shown that there exists a unique term g(x a ,p 0 ), 

, . 21 (mi + to 2 )toiTO2 , 

,(x«,P») = ¥l {Wl)lrh ( n 12 Pi) 

1 


r 12 

4 


d -3 


2. Integral of h 2 ^ 


0{d — 3). 


( 8 . 10 ) 


The integral of the density hd 2 ) is not convergent at 
spatial infinity and it also develops UV poles. After mak- 


/ 97 m 2 

2 n 

\ 20mi 

16/ 
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+ 


+ 


/ 217711 

V^6“ 



687m 2 \ 


40 


j (n 12 • Pi) (n 12 • p 2 ) 


237?7l 2 ^ 


40mi 


77i 2 (n 1 2 • Pl)Pi 


+ 

+ 

+ 


7777! 

"48" 


+ 


/ 15977772 

v 120 


307777,2 ^ 
40 J 

77771 

— 


(n 12 • P2)P? 
)(ni2-Pi)(Pi P2) 


(1 -o- 2), 


( 8 . 12 ) 


such that 


xli 2 ) 4 (*a,Pa) + -^g(x a , Pa) = F{x a ,p a ), 


d t 


(8.13) 


where 


F(x a ,Pa) 


2 M r . T . , 2 
5(167t) 2 ^ ^ 


(8.14) 


Here is the Newtonian quadrupole moment of the bi¬ 
nary system, 


Therefore as the contribution of the density /7d 2 ) to the 
4PN Hamiltonian we take 

Hlp e N 2 ' 2 (x a ,p a ; d; s, C ) = x 2 { i 2 )(^a,Pa) 2 2 3 

. 2,2,2/ r 12 X(12) ( X “>Pa) 

+ X ( 12 )(x a ,Pa)ln ir + d 3 

+ F(x a , p a ) (In ^ + C) + 0(d - 3). 

(8.19) 


Let us finally mention that the integral of the last term of 
the density /i 2 { 2) , B^jA^hJ^, though IR divergent, 
can be in fact regularized uniquely. This is so because 
the difference between results of different methods of reg¬ 
ularizing this term is a total time derivative. 


3. Integral of h^f 2 ) 


The density hf : \ 2 ) an exac t divergence, it can be writ¬ 
ten in the form 

KdSbSb * s 5)«< A ; 1 '>S«>) = J**. &*» 

where 


hj .= ^2m a (x z a x J a - (8-15) 

a 

All time derivatives on the right-hand side of the formula 
(18.1411 were eliminated by means of Newtonian equations 
of motion. By virtue of Eq. (18.131) the result of Eq. (18.1111 
we can rewrite in the following form: 


E k := C^{A- d 'hJl A „) - C^ k (A~ d v hJ^). (8.21) 

The surface integral associated with E k is not convergent 
at spatial infinity in d = 3 dimensions. However in d 
dimensions E k ~ r 5_2d for r 00 , so the surface integral 
behaves like r 4_d and it vanishes in the limit r —> 00 
when d is large enough. Therefore this term does not 
contribute to the Hamiltonian. 


/ 


d d x 


i2, 2 

n {12) 


+ 


d t 


(l(Xa,Pa) 


In 


ri2 \ 


2,2,1 

X(l 2 ) a mb 


( x oj Pa; 


C) 


D. Removing UV poles 


2,2,3/ \ 

T 12 X(12) ( x a;Pa) 


1 2,2,2/ \| ' i/i 1 

+ X(i 2 ) ( X a; Pa) hi —= T 


d -3 


ri2 


where 
2,2,1 


+ F(x a , p a ) In — + 0(d — 3), 
s 


(8.16) 


^(12)amb( Xa ’ P fl ’ ^(12)amb( Xa ’ P fl ’ 

d 


g( x a,Pa)-^(ln^y (8.17) 


We have found (see Appendix IA 311 that the ambiguity 
of the term X(i 2 )amb( Xa ’ Pa! C) can be expressed, up to 
adding a total time derivative, as a multiple of the term 
F introduced in Eq. (18. 141) . 

X(12)amb( X °’P a ’ = X(l 2 ) ( x a;Pa) + C-F(x 0 , p a ) 

+ (total time derivative). (8.18) 


As the result of UV and IR regularizations described 
in the previous subsections, we have obtained the 4PN 
near-zone Hamiltonian of the following form: 

7f4P S N (x a , p a ; d; S, C) = -ffJpN^Xa, p a ; d) 

+ #Ip S N 2 , 1 ( x a, Pa; d) T 774p g N 2 ’ 2 (x a , p a ; d; s, C) 

- / \ - f M ? ’12 . X3( x a,Pa) 

= Xl( x a, Pa) + X2 ( x a ; Pa ) hi -j- + -- 

d — 3 

+ E(x a ,p a )(ln^ + C) + 0(d-3), (8.22) 

where 

Xfe( x a,Pa) := X(( 2 )( x a,Pa) +X(12) (XoiPo) 

+ X(i 2 ) fe ( x a, Pa), k = 1, 2, 3. (8.23) 

In the next step we remove both the DR-related scale £q 
and the pole terms proportional to 1/ (d — 3) by adding a 
total time derivative. To do this we have found a unique 
function Z?(x Q ,p a ;d), 
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L>(x a ,p 0 ;d) = 


m i 


(16tt ) 3 rf 2 


( X 39 2 \ 

( [( n i2 • Pi ) 2 - Pi] ( n i2 • P2) + 2 (ni 2 • Pi)(pi • P2)} + g(ni2 • p 2 )p 2 j 


d -3 


+ [pi - ( n i 2 ■ P 1 ) 2 ] (ni 2 ■ p 2 ) - 2(ni2 • Pi)(pi • p 2 )} - 2(ni 2 ■ p 2 )p^ In 3 j + (le 2), 


(8.24) 


such that the sum 

#4P S N (x a ,p a ;d;s,C) + ^L>(x a ,p a ;d) 

has the finite limit as d —> 3. This limit we take as the 
regularized value of the 4PN near-zone Hamiltonian: 

ff near- zon e (s) ( X(j ; p a; Q) : = lim (#][p S N ( x a, Pai d] S, C) 

d —>3 V 

+ ^ £) ( xq ’ p a ’ d ))' ( 8 ‘ 25 ) 

This Hamiltonian depends on the IR regularization scale 
s and on the dimensionless constant C and can be rewrit¬ 
ten as 


+ ^ sym(s) [x a , Po ], (8.28) 

where we have used brackets [•,•] to emphasize that the 
tail Hamiltonian H^p^ synl ^ is a functional of phase- 
space trajectories x a (t), p a (t). Reference [^3| also com¬ 
puted the value of the constant C, 


C = — 


1681 

1536’ 


(8.29) 


and showed that the total 4PN Hamiltonian (18.2811 does 
not depend on the scale s and can be written as 


i? 4 PN[x a , p a ] — ^ipN ( x O) Pa) + -f^4PN° Cal [ x a ) Pa]) 


(8.30) 


where the local piece of the 4PN Hamiltonian reads 


„near-zone («), . — fj local 0/ \ 

rI 4PN l x O) Pa> ) — n 4P N t x a,PaJ 

+ F(x a ,p a )(ln^ + c'j, 

(8.26) 

where the uniquely computed part H^§ 1 0 of the near¬ 
zone Hamiltonian reads 

# 4 PN 10 ( X a)Pa) == Xl ( x a, Pa) + Hm f %2 (x a , p a ) In ^ 

d -¥3 y to 

, X3(x a ,Pa) , d \ 

+ d -3 +t t D <X:P‘^)- (8-27) 

Let us stress again that H\ pg 1 0 is without DR-related 
scale and poles in 1 /(cZ — 3). 


H \PN ( X a, Pa) = ^4PN °( x a, Pa) + CF(x a , p a ), (8.31) 

and the nonlocal-in-time piece can be written as (from 
now on we restore the constants c and G) 

rrnonlocalp. ] _ 1 G M ... 

- n 4PN l x a: PaJ — ~ g — fs 1 t? 

r +°°... du 

X Pf 2ri2 / C / Iij(t + V)y~. (8.32) 

Here 1 denotes a third time derivative of the Newto¬ 
nian quadrupole moment of the binary [defined in Eq. 
(18.151) ] and Pf 2ri2 / C is a Hadamard partie finie with time 
scale 2?'i 2 /c [see Eq. (4.2) in hi] for the definition of the 
Pf operation]. 

The total 4PN-accurate conservative matter Hamilto¬ 
nian is the sum 


E. Total 4PN-accurate conservative matter 
Hamiltonian 

Reference 51] showed that the total 4PN conserva¬ 
tive matter Hamiltonian is the sum of the local-in-time 
near-zone Hamiltonian (18.2611 and the time-symmetric 
but nonlocal-in-time tail Hamiltonian ^, 

tt r 1 Tjnear-zone (s) / 

#4PNl x a,Pa] = # 4 p N '(x a ,p a ;C7) 


#<4PN[ x a,Pa] = ^ m a (? T H N (x a , p a ) 
a 

+ i?l PN (x a , Pa) + iJ 2 p N ( x a, Pa) 

T Lf 3 pN ( x a, Pa) T L/ 4 PN [ X a, Pa] • (8.33) 

The explicit formulas for the local-in-time Hamiltonians 
from the Newtonian up to the 4PN level [the local piece 
(18.311) of the 4PN Hamiltonian given below incorporates 
the value (18.291) of the constant C], valid in the generic, 
i.e., noncenter-of-mass, reference frame, are as follows: 
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2 mi 2 ri2 


(8.34a) 


c H 1PN (x a ,p a ) = - -—+ -- 

8 mf 8 r 12 


1 (P?) 2 , 1 Gm 1 m 2 f p? , 1 ,(pi-p 2 ) , 0 (n 12 • pi)(n 12 • p 2 ) 


-12 + 14 ^ + 2 ■ 

mf m.\m .2 


1 Gmim 2 G(mi + m 2 ) 
4 r 12 r 12 


+ (l O 2), 


(8.34b) 


C 4 l?2PN(x a , p a ) = 


1 (pf) 3 lGm 1 m 2 ( c (P?) 2 HPiPl (Pi ’ P2) 2 , _ Pi (ni2 • P2) 2 


16 m\ ' 8 ri 2 \ m\ 2 m\m\ m 2 m| 


(pi • p 2 ) (n 12 • pi)(ni 2 • p 2 ) 3 (n 12 • pi) 2 (ni 2 • p 2 ) 5 


1 G 2 mim 2 

4 r 12 


m 2 (1„4 + l4) _ V + m2) 27(P.-K) + 6(.ii 2 Pi)(n, iK ) 

\ mf mf J 2 mi m2 


1 Gffliffl 2 G 2 (mf + 5miTO 2 + m 2 ) 
8 r 12 r? 2 


+ (l^ 2), 


(8.34c) 


r „ _ n 5 (P?) 4 , 1 Gmim 2 ( (p?) 3 . ((Pi ' P2) 2 + 4p 2 p 2 )p 2 

C ^3PN(Xa, P a)-- — + -14—^+4-- 


Pi ( n i2 • Pi) 2 (ni2 • p 2 ) 2 (pf (ni 2 ■ p 2 ) 2 + p 2 (ni2 • Pi) 2 )Pi p 2 (pi • p 2 )(ni 2 • Pi)(ni 2 ■ p 2 ) 


P? (Pi ' P2)(ni2 • P2) 2 , (”p?p 2 - 10(pi • p 2 ) 2 )(ni 2 • pi)(ni 2 • p 2 ) (p? p| - 2 (pi • p 2 ) 2 ) (pi • p 2 ) 


, , c (Pi • P 2 )(ni 2 • Pi) 2 (ni 2 • p 2 ) 2 1 o pf (n 12 • pi)(ni 2 • p 2 ) 3 | c (ni 2 • pi) 3 (ni 2 • p 2 ) 3 

' 10 1 1 13 (( *3 „,3 „,3 ' 0 ^ 3^,3 


G 2 m 1 m 2 ( 1 0 -7 „„n(Pi ) 2 115 „„ Pi (Pi • P 2) , 1_ 25 (pi • p 2 ) 2 + 371 pf pf 

2 I 1 a (^1 27777-2) 4 „ mi 2 + - q 7772 2 9 

r^ 2 \ 16 m\ 16 mjm2 48 mfmf 


17 p 2 (n 12 • pi) 2 5 (n 12 ■ pi) 4 _ 1^ (l5pf (n 12 ■ p 2 ) + 11 (pi • p 2 ) (n i2 • Pi))(n 12 • pi) 

16 m 3 12 to 3 8 m\m 2 


3 (ni 2 • pi) 3 (ni 2 • p 2 ) 125 (pi • p 2 ) (ni 2 • pi)(ni 2 • p 2 ) 10 (ni 2 • pi) 2 (n 12 • p 2 ) 2 

777-1 3 + i 0 ^2 2 2 o ^2 2 2 

2 m{m2 12 mfmf 6 mfmf 


- 1(220™, + 193m 2 ) p|(n 'rf )2 ) + - 1 

48 mfmf / r\ 2 l 48 


^425 to 2 + ^473 — ^ 71 ‘ 2 ) mim 2 + 150 m 2 ^ 


fn’j( 2 , 2 ( , L, 1 2 \ \ (Pi ' P2) . 1 / on 2 , 3 2 \ \ (n 12 • p 

■( 77(mj+m 2 ) + 143--7T mim 2 )-I- — ( 20 m, - 43+ -tt mim 2 I-5- 

1 \ “V 4 / J m.im.2 16 \ V 4 / J mf 


+ — ( 21(mj + mf) + (119 + 


2 )mim 2 ^ 


(n 12 • pi)(ni 2 ■ p 2 ) \ lG 4 mim|//227 21 9 \ 

- +fi-4- “5-T 7r TO 1+ TO 2 

mim 2 / ® r i2 \ V 3 4 ) 
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+ (l^ 2), 


(8.34d) 


rrlocal ( 
il 4PN l 

X-a 1 Pa) — 

_l_ 

G 3 mirri2 

1 

r 3 
' 12 

_L 

G 4 ? 7 li? 7 l 2 

> 4 

r 12 


7(Pi) 5 , Gtoi? 71 2 , , G 2 mim 2 , , 

956m 9 H - — - #4 8 (x 0 ,p a ) H -- TOi H i 6 (x a , p Q ) 

(m? 7J 4 4i(x a , p a ) + TO 1 TO 2 #442 (x a , Pa)) 

(to.? #42l(x a ,Pa) + TO?TO 2 #422(x a , Pa)) + G ™ 1??72 g 40 (x a , p a ) + (lH 2), 

\ / Tin 


(8.34e) 


„ , 45(p?) 4 9(ni2 • pi) 2 (ni 2 ■ p 2 ) 2 (p?) 2 15(ni 2 • p 2 ) 2 (p?) 3 9(ni 2 ■ Pi)(ni 2 ■ p 2 )(p?) 2 (pi • P 2 ) 

48 “>P° 128m? 64 to?to| 64m?m 2 16 to?to 2 

_ 3(p?) 2 (pi ■ p 2 ) 2 + 15(ni2 ■ pi) 2 (p?) 2 p? _ 21(p 2 ) 3 p 2 _ 35(ni 2 ■ pi) 5 (ni 2 ■ p 2 ) 3 
32 toiTO? 64to.iTO? 64toiTO? 256?7i?TO2 

25(ni2 • Pi) 3 (n 12 ■ p 2 ) 3 p? 33(ni 2 • pi)(ni 2 • P 2 ) 3 (p?) 2 _ 85(ni 2 ■ pi) 4 (ni 2 • P 2 ) 2 (Pi • P 2 ) 

128to???t.2 256toi?tt.2 256to?TO2 

_ 45(ni2 • Pi) 2 (n 12 ■ p 2 ) 2 p?(pi • P 2 ) _ (n 12 • P 2 ) 2 (p?) 2 (Pi • P 2 ) 25(ni 2 • pi) 3 (ni 2 • p 2 )(pi • P 2) 2 

128to?to 2 256to?to| 64m?TO| 
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64 to?to| 64to?to? 64TO 4 m? 
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256 m\m\ 128m?m 2 256m?m 2 
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64 m|m 2 64mfm 2 4infrri2 

(n 12 • P 2 ) 2 P?(Pi • P 2) 2 _ 5(n 12 • pi) 4 (ni 2 • p 2 ) 2 pj 21(n 12 • Pi) 2 (n i2 • P 2 ) 2 P?pj 
16 to?to. 2 64rn 4 TO 2 64 rn 4 m 2 

3(n 12 • p 2 ) 2 (p?) 2 pl (n i2 • Pi) 3 (ni2 • p 2 )(pi • P2)P? , (n 12 • P!)(n 12 • p 2 )p?(pi • p 2 )p? 


32 mfmf 


4 TO 4 m 2 
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16m?m 2 32mfm 2 64 to?to. 2 

3(n 12 -pi) 2 p 2 (p 2 ) 2 7(p 2 ) 2 (p 2 ) 2 

32 mfmf 128m 4 m 2 


32 m 4 mf 


(8.34f) 


#46(x a ,Pa) = 


369(n 12 • pi) 6 889(n 12 • Pi) 4 p 2 , 49(ni 2 • pi) 2 (p?) 2 63(p?) 3 


160m? 


192m? 
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1280to 4 to 2 480to 4 to 2 3840to 4 to 2 
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I 6 OO 777 2 
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^422(x a ,Pa) 
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(8.341) 


For completeness let us supplement the generic 4PN- 
accurate local Hamiltonian given above by its center-of- 
mass expression. The center-of-mass reference frame is 
defined through the condition 


Pi + P2 = 0. (8.35) 


Let us introduce the reduced mass /i of the system and 
its symmetric mass ratio u, 


TOiTO 2 /i TOlTO 2 

M ' M (mi + m^) 2 


(8.36) 


which is the sum of different PN contributions, 


H< 4pn[i", p] = H N (r, p) + H 1 PN (r, p) + H 2 P n( r, p) 

+ ^3P N (r,p) +F7 4 PN[r, p], (8.39) 


where the 4PN Hamiltonian is the sum of local- and 
nonlocal-in-time parts, 


It is convenient to introduce the reduced variables 
xi 2 Pi £2 


r := 


GM' 


Pi 

p := — = 


(8.37) 


#4PN [r, p] = #4 pn (r, p) + 77 4 n pN° ca [r, P] (8-40) 


(with r := |r| and n := r/r). We also define the reduced 
4PN-accurate center-of-mass Hamiltonian 

H< 4PN [r, p] := g ~ 4PN[r ’ P] ~ Mc2 , (8.38) 

[X 

_I 


The local Hamiltonians from to H\ p^ 1 are equal to 
(let us recall that their coefficients depend on masses m\, 
m 2 only through the symmetric mass ratio v) 


Hn (r,p) = - -, (8.41a) 

2 r 

C 2 H 1PN (r, p) = i(3^ - 1 )(p 2 ) 2 - i|(3 + p)p 2 +p(n • p) 2 }i + ^- 5 , (8.41b) 

c 4 ^ 2 PN (r, p) = ^ (1- - + 5p 2 ) (p 2 ) 3 + (5 - 20^ - 3 p 2 ) (p 2 ) 2 - 2p 2 (n • p) 2 p 2 - 3i' 2 (n • p) 4 }^ 

+ i{(5 + 8p)p 2 + 3p(n ' P) 2 }^ - j(l + (8.41c) 

c 6 H 3PN (r, p) = 2^ (-5 + 35^ - 70^ 2 + 35 p 3 ) (p 2 ) 4 + ^ 7 1 (-7 + 42i2 - 53 p 2 - 5p 3 ) (p 2 ) 3 
+ (2 - 3p)i/ 2 (n • p) 2 (p 2 ) 2 + 3(1 - v)v 2 [n ■ p) 4 p 2 - 5p 3 (n • p) 6 ji 

+ {^ (—27 + 136p + 109p 2 ) (p 2 ) 2 + ^(17 + 30i/)i/(n • p) 2 p 2 + -L(5 + A3u)u(n • p) 4 |l 
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(8.41e) 


IX. THE POINCARE INVARIANCE 


There are several possibilities for partial checks of our 
results: first the test-body limit and second the linear in 
G part of the 4PN Hamiltonian through comparison with 
the post-Minkowskian results achieved in Ref. [54]. Both 
these checks were already performed by us in Ref. [52]. 
The most important tool however is the Poincare invari¬ 
ance, discussed for the first time in the context of the 
ADM Hamiltonian approach to the two-body problem 


in Ref. [35[ . Poincare symmetry holds because our (iso¬ 
lated and conservative) two-point-mass system is living in 
asymptotically flat spacetime with its globally conserved 
quantities: energy H , linear momentum P, angular mo¬ 
mentum J, and Lorentz boost vector K = —P t + G, 
where G denotes the center-of-(mass)energy vector. All 
these quantities are realized as functions on the two-body 
phase space (xi,x 2 ,pi,p 2 ), whose usual Poisson brack- 
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ets, 


{/(x 0 ,p a ), 5 (x a ,p a )} := 

a 


f df dg 
V dx l a dpai 


df dg\ 
dpa i dxi ) 


(9.1) 

satisfy the Poincare algebra relations, 

{Pi ,Pj} = 0, { Ji , Jj} = CijkJk , (9.2a) 

{Ji, Pj} = eijkPk, (9.2b) 

{Pi,H} = 0, {J ly H} = 0, (9.2c) 

{Ji, Gj } = CijkGk, (9.2d) 

{G h H}=P t , (9.2e) 

{( G l ,P 0 }=c~ 2 H5 lj , (9.2f) 

{Gi, Gj} = -c- 2 e ijk J k . (9.2g) 


These relations have to be fulfilled through 4PN order. 

The total linear and angular momenta have universal 
forms, 


f)(^a, Pa) — ^ figi ■ <l«(x a ,p a ) — ^ ] ^ijk^gPak, 
a a 

(9.3) 

and they exactly satisfy Eqs. (I9.2al) (I9.2cl) . We will con¬ 
struct the boost vector G as a three-vector from x a and 
p 0 only; therefore, the relation (I9.2dl) will also be ex¬ 
actly satisfied. Consequently Eqs. (I9.2cl) ( |9)2gl ) are the 
only nontrivial relations which have to be satisfied by the 
vector G. The Hamitonian H entering Poincare algebra 
(EH) is the full 4PN-accurate Hamiltonian, 

#<4PN[x a , Pa] = ^<4 Pn( X 0) Pa) + #4PN° Cal [ x a, Pa], 

(9.4) 

where the local-in-time part reads 


^< 4 Pn( X o, Pa) — 'y ] iTlgC 2 + Hj<i(x a , p a ) 

a 

+ 7Ji PN (x a , p a ) + Ef 2 PN(x a , Pa) 

+ ^?3Pn(x 0 , Pa) + l?4p C N ( x a. Pa)- 

(9.5) 


Because the nonlocal-in-time piece -Hipn 008,1 is Galileo 
invariant [see Eq. (5.15) in j5lj for the proof], it is enough 

I 


to restrict the 4PN-accurate Hamiltonian to its local part 
H < 4 Pn when looking for the 4PN-accurate boost vector 

G7 

We have found the 4PN-accurate boost vector G using 
the method of undetermined coefficients employed at the 
3PN level in Ref. [35]. The generic form of the three- 
vector G reads 

G (x a , Pa) = Y ( M 0 (xt,,p f) )Xa+fV Q (x & ,p & )pa), (9.6) 


where the scalars M a and N a possess the following 4PN- 
accurate expansions 


M a =m a + M a 1PN + M 2PN + M 3PN 
N a = N 2PN + 1V 3PN + N* pn . 


M. 


4PN 


(9.7a) 

(9.7b) 


Let us note that M a and N a reduce to m a and 0, re¬ 
spectively, in the Newtonian approximation. Next we 
write the most general expressions for the successive PN 
approximations to the functions M a and N a as sums of 
scalar monomials of the form 


Cn 


r 


—no 

12 



( Pi • P2 

V rnirri2 


n 2 



n 3 


n !2 • Pi 

mi 


n i2 P 2 
m 2 


m, u m 9 


(9.8) 


where no,..., nr are non-negative integers. To con¬ 
strain the possible values of no,... ,nr we employ di¬ 
mensional analysis, Euclidean covariance (including par¬ 
ity symmetry), and time reversal symmetry (which im¬ 
poses that M a is even and N a is odd under the operation 
Pa —» —Pa)- We also use the 1 O 2 relabeling symme¬ 
try. At the 4PN level the most general pattern for the 
functions M 4PN and ./V 4PN involves 210 dimensionless 
coefficients c n . 

To find the functions M 4PN , ..., M 4PN and 7V 2PN , 
..., 1V^ PN it is enough to use Eq. (I9.2dl only. The 3PN- 
accurate functions were constructed in Ref. [35[ (for com¬ 
pleteness we give below their explicit expressions). At 
the 4PN level the relation (I9.2cll yields 525 equations to 
be satisfied. We have found a unique solution to these 
equations [and we have then checked that this solution 
satisfies the remaining Poincare algebra relations (I9.2fl) 
and Q9.2g|) ]. The explicit forms of the functions M 4PN , 
..., M 4P ” anc j jv 2PN , ..., 7V a 4PN read 


1 pj 1 Gm\m2 


C 2 M 1 1PN (x a , Pa ) = ^~ 
C 4 M 1 2PN (x a ,p a ) = - 


2 mi 2 n 2 
1 (p 2 ) 2 1 GmiW2 ( 


r 12 


(9.9a) 


8 mf 4 ri 2 \ 

1 Gm\m 2 G{m\ -(- m 2 ) 


- 5 


Pi_ _ P| 7 (Pi ' P 2 ) (ni 2 • Pi)(ni 2 • £ 2 ) 

771? 777? 777 x 7772 777 1 777-2 


4 r 12 


(9.9b) 
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c 6 M^(x a ,p a ) = -^f + - 


1 (pf ) 3 1 GTO1TO2 ( (Pi ) 2 , (pi ) 2 -- p?pi (Pi ' P 2 ) 2 , o Pi ( n 12 • P2) 


16 m\ 16 ri 2 \ m 4 m: 


— 1 1 J - _ o Ifli 1 - 

-Li 9 9 " 9 9 

mi mi mi mi 


pi (ni 2 • Pi ) 2 (pi • p 2 ) (ni 2 • Pi)(ni 2 • p 2 ) 0 (n i2 • Pi) 2 (ni 2 • p 2 ) 2 

I 1 22 ^" 22 ^ 2 2 


1 G 2 ?77 i? 77 2 

+ 24 W 


— (31mi + 5m 2 ) 


(112toi + 45to 2 )-^4 + (15rni + 2 m 2 )-J-4--(209mi + 115m 2 )-^-^—^ 

mf 777-2 2 777-lTO 2 

(ni 2 • Pi)(n i2 • p 2 ) (n i2 ■ pi) 2 _ (n i2 • P 2) 2 \ _ 1 Gmim 2 G 2 (m\ + 5mim 2 + mi) 
mim 2 mi 777 2 / 8 r i2 r ? 2 


c 8 M 1 4PN (x a ,p a ) = - (x a , p a ) + G ” 2 1??l2 (mi Af 44 i(x a , p a ) + m 2 M 442 (x a , p a )) 


G 2 777 -i 777. 2 


128mi ri 2 


G 3 mim 2 


(m 2 M 4 2 i(x a , p a ) + mim 2 M 422 (x a , p a ) + mi M 423 (x a , p a )) 


G 4 mim 2 ,, . . 

—-j-M 40 (x a ,p a ), (9.9d) 


, \ _ 13 (pj ) 3 _ 15 (n 12 ■ pi) 4 (ni 2 ■ p 2 ) 2 45 (ni 2 ■ pi) 2 (n 12 • p 2 ) 2 pj _ 91 (ni 2 ■ p 2 ) 2 (pj ) 2 

46 °’P a 32777 ® 256 ? 7 ifmi 128 m 4 ? 7 ii 256m\rn^ 
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32 m 4 777 -i 32 m 4 mi 64 mfmi 
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32 m 3 ?772 32 m 3 ?772 32 m 3 777 3 
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_ 7(ni2 ■ Pi)(ni2 • p 2 )p| _ 7(pi ■ p 2 )p| _ 9(p^) 2 
6mim 2 48mim 3 32m 2 ’ 
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(9.9h) 


c 4 7v 2PN ( XaiPa ) = —- G (n 12 -p 2 ), 

c 6 lV 1 3PN (x 0 ,p a ) = ^ ° (2 (pi • p 2 )(n 12 • p 2 ) - p 2 (n 12 • Pi) + 3(n 12 • pi)(n 12 • p 2 ) 2 ) 

8 777-1777-2 V / 
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The fulfillment of the Poincare algebra does not imply 
a complete check of the Hamiltonian, but rather a check 
of all terms besides the purely static ones. Of course, the 
Poincare algebra is invariant against canonical transfor¬ 
mations, particularly those induced by coordinate-gauge 
transformations, so the single terms in the Hamiltonian 
have no direct physical meaning. The reader interested 
in a representation of a higher-order PN Hamiltonian 
through center-of-mass and relative coordinates is re¬ 
ferred to |72l|. 
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Appendix A: Regularization 

In this appendix we describe techniques which we have 
used to regularize divergent integrals which appear in our 
paper. 

1. Three-dimensional Riesz-implemented 
Hadamard regularization 

The Riesz-implemented Hadamard (RH) regulariza¬ 
tion was developed in the context of deriving PN equa¬ 
tions of motion of binary systems in Refs. [32|, [4l[ (see 
also 0 ) to deal with locally divergent integrals com¬ 
puted in 3 dimensions. The RH regularization relies on 
multiplying the full integrand, say *(x), of the divergent 
integral by a regularization factor, 
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and studying the double limit ei —> 0 , e 2 —> 0 (here Si 
and s 2 are arbitrary three-dimensional UV regularization 
scales). Let us thus consider such an integral performed 
over the whole space R 3 and let us assume that it devel¬ 
ops only local poles (so it is convergent at spatial infin¬ 
ity). The value of the integral, after performing the RH 
regularization in 3 dimensions, has the structure 

/ RH (3;e 1 ,e 2 ) := f i(x)^ ^ d 3 z 

Jr3 v si / V s 2 / 


(1 

1 1*12 > 

( 1 

, r 12 \ 

— 

+ In- 

+ C2 — 

+ In — 

Vei 

Sl ) 

Ve2 

S2 ’ 


+ 0(e i,e 2 ). (A2) 

In the case of an integral over R 3 developing poles only 
at spatial infinity (so it is locally integrable) it would be 
enough to use a regularization factor of the form (r/ro) e 
(where ro is an IR regularization scale), but it is more 
convenient to use the factor 


= A~, 


(a + b)e 


In 


r 12 
ro 


0(e). 

(A4) 


Many integrals appearing in Eqs. (IA2I) and (IA4I) we 
compute using three-dimensional form of the following 
formula, first derived in d dimensions probably by Riesz 
0: 



0 

2 


d d x = 7 r d / 2 


a+(3+d 

' 12 


r(^)r( 

k 2 > 

in 


r(-f)T 

H) 

1 r ( 

f a+0+2d^ 


(A5) 


To compute the 4PN-accurate Hamiltonian one needs 
to employ a generalization of the three-dimensional ver¬ 
sion of the Riesz formula (IA5I) for integrands of the form 
r“ r% s 7 , where 


(-) 0£ (-) 6e (A3) 

Vr 0 / Vr 0 / 

and study the limit e —> 0. Let us denote the integrand 
again by i(x). The value of the integral, after performing 
the RH regularization in 3 dimensions, has the structure 

/ RH (3;a,&,e):= [ i(x) (^)‘“7^)* d 3 z 
J R3 V ro / v r 0 / 


s := ri + r 2 + r 12 . (A 6 ) 

Such formula was derived in Ref. [H] and it reads 

J r“ s 7 d 3 a; = R(a, f3, 7 ) r“ 2 f/3+7+3 ) (A7a) 

where [let us note that the formula given below is invari¬ 
ant under the permutation (a -o- b)] 


R(ot, (3, 7 ) := 2t r 


T (a + 2 ) T (/3 + 2 ) T (-a- 

r (—7) 


/3 - 7 - 4) 


x [J 1/2 (a + 2, -a - 7 - 2) + J 1/2 (/3 + 2, -/3 - 7 - 2) - / 1/2 (a + /3 + 4, -a - /3 - 7 - 4) - l] . (A7b) 


The function /i / 2 in Eq. (IA7bl) is defined as follows: 


h /2 (x,y) 


B 1/2 (x,y) 
B (x, y) 


(AS) 


where B stands for the beta function (Euler’s integral of 
the first kind) and B 1 / 2 is the incomplete beta function; 
it can be expressed in terms of the Gauss hypergeometric 
function 2 F\\ 


B 1/2 (x, y) = ^ 2^1 - y, x; x + 1; 0 . (A9) 

The regularization procedure based on the formulas 
(IA5D and (IA71) consists of several steps. We enumerate 
them now. The most general integrand we have to con¬ 
sider has the form 


(rii • Pi ) 91 (n 2 • pi ) 92 (ni • p 2 ) 93 (n 2 • p 2 ) 94 ri r^s 7 , (A10) 


where q ±, ..., q± are non-negative integers and 7 is a 
negative integer. We first eliminate the unit vector n 2 
by the identity 

ri , r !2 / A11 i 

n 2 = —ni -|-n 12 . (All) 

r2 r 2 

We thus plug Eq. (IA11D into (IA10I) and expand the scalar- 
product ni • ni 2 by means of the relation 


ill • ni 2 


r 2 _ 2 _ 2 
'2 '1 '12 

2riri 2 


(A12) 


After this the most general integrand reduces to 


(ni • Pi) 9l (ni • p 2 ) 92 rirfs 7 , (A13) 


where again qi and q 2 are non-negative integers. 

We perform integration in prolate spheroidal coordi¬ 
nates. By using these coordinates it is possible to rep¬ 
resent integrand (IA13I) as a linear combination of inte¬ 
grands of the type rJV^s 7 . To show this let us locate the 
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particles in the focal points of the prolate spheroidal co¬ 
ordinates [they lie along the z axis of the Cartesian coor¬ 
dinate system (x,y,z)\, so the particles’ position vectors 
have the following Cartesian coordinates: 

Xl = (0,0, -r 12 /2), x 2 = (0,0, r 12 /2). (A14) 

The Cartesian components of the unit vector ni := (x — 
xi) /r± can be expressed by t \, r 2 and the azimuthal angle 


where ( B) is the average over the angle <fi, 

{B):= hC (A19) 

After this step the integrand (IA13I) becomes the linear 
combination of the type 

(A20) 



COS0, 

(Al5a) 


sin 0, 
(A15b) 


n 


Z 

1 


r 2 i r 2 _ 2 

; 12 T Tj_ 1 2 

2 r 1 r 12 


(Al5c) 


Without loss of generality we can place the vector pi 
in the (x,z) plane (we can also assume that p\ x > 0). 
One can then show that the Cartesian components of 
the vector pi are as follows: 


Pix = y^P? - (ni 2 • Pi) 2 , (Al6a) 

Piy = 0, (A16b) 

Pn = —(n 12 ■ Pi). (A16c) 


The x and z Cartesian components of the vector p 2 read 


P2x = 


(Pi • P 2 ) - (n 12 ■ Pi)(ni 2 • p 2 ) 
\/Pi - ( n i 2 • Pi) 2 


(Al7a) 


P2z = — (ni 2 • p 2 ); 


(A17b) 


the y component of the vector p 2 can be computed from 


P2y = ± 


Pi 


Plx 


■pIz 


(A17c) 


We have checked that the result of the procedure de¬ 
scribed below does not depend on the ± ambiguity in 
Eq. (IA17cl) . 

Making use of Eqs. (IA15D (IA17I) we compute scalar- 
products (ni • pi) and (ni • p 2 ) in Eq. (IA13I) . After this 
the integrand (IA13I) becomes a sum of terms of the form 
A(ri, r 2 )-B(</>), where R(</>) is a polynomial in sin</> and 
cos </>. Each of these terms we integrate in the following 
way: 


A(n,r 2 )B((j))d 3 x = / A(n,r 2 ) d 2 x / !?(</>) d<^> 


(A18) 


where the constant coefficients c/ may depend only on 
?T 2 , P?, (P i • P 2 ), p|, (n 12 • pi), and (n 12 • p 2 ). The 
integral of (IA20D is computed by means of Eq. (IA2II or 
(El and with the usage of formulas El and El- 
Appendix A6 of Ref. contains generalization of the 
above presented procedure to d space dimensions (and it 
employs prolate spheroidal coordinates in d dimensions). 


2. UV corrections 

Reference [32] showed that the three-dimensional RH 
regularization described above used to derive the 3PN 
two-point-mass Hamiltonian gave ambiguous results. 
Namely, by means of integration by parts (assuming 
that all involved integrals are convergent at infinity, so 
all surface terms can be neglected) one can replace one 
form of Hamiltonian density (or its part) by some other 
form. Integration of both equivalent densities should 
give the same result, but it did not. To correct the re¬ 
sult of the three-dimensional RH regularization (i.e., to 
remove ambiguity), Ref. (3l] (see Secs. 3 and 4 there) 
developed dimensional-regularization (DR) technique^ 
which we have also used to make the results of the three- 
dimensional RH regularization of the locally divergent 
part of the 4PN-accurate Hamiltonian unique. 

The technique of Ref. [3l| boils down to the computa¬ 
tion of the difference 

lim H l ^(d) - i7 4 R p H N loc (3), (A21) 

a—>-3 

where i7|^ loc (3) is the “local part” of the Hamilto¬ 
nian obtained by means of the three-dimensional RH 
regularization [it is the sum of all integrals of the type 
J RH (3;ei,e 2 ) introduced in Eq. (|A2[) [: H l ^ N (d) is its d- 
dimensional counterpart. 

Reference |3l| showed that to find the DR correction to 
the integral / RH (3;ei,e 2 ) related with the local pole at, 
say, x = xi, it is enough to consider only this part of the 
integrand *(x) which develops logarithmic singularities, 
i.e. which locally behaves like 1/rf, 

i(x) = • • • + Ci(ni) rf 3 + • • • , when x —» xi. (A22) 


3 The presentation of this technique given below is an improved 
and more complete version of the explanations contained in Sec. 
Ill of m. 


= (B) / A(ri,r 2 ) d 3 x, 


































33 


Then the pole part of the integral (IA2I) (related with the 
singularity at x = Xi) we recover by RH regularization of 
the integral of ci(ni) rf * * 3 over the ball D(xi,t'i) of radius 
i\ surrounding the particle xi. The RH regularized value 
of this integral reads 



where C\ is the angle-averaged value of the coefficient 
ci(ni). The expansion of the integral d RH (3;ei) around 
ei = 0 equals 

/^ H (3; ei) = ci ("-b In — \ 0(e±). (A24) 

Vei s\J 

The idea of the technique developed in ( 3l] relies 
on replacing the RH-regularized value of the three- 
dimensional integral J RH (3;ei) by the value of its d- 
dimensional version Ii(d). We thus consider the d- 
dimensional counterpart of the expansion (IA22I) . It reads 

i{x) = -b^o <d 3 ^i(d; ni) rf~ 3d H-, when x —» Xi, 

(A25) 

where is the scale which relates the Newtonian Gn 
and the D-dimensional (D = d + 1) Go gravitational 
constants, 

G D = G N £ d ~ 3 * * * * * . (A26) 

The number k in the exponent of £^ d 3 "* is related with 
the momentum-order of the considered term [the term 
with k is of the order of O(p 10 ~ 2k ), where k = 1,..., 5]. 
The integral Ii(d) we define as 

h (d) := e k 0 {d ~ 3) [ Ci (d; m) rl~ 3d d d n 

J B(x i,ti) 

= tf (d - 3) a(d) r^dn, (A27) 

Jo 

where Ci (d) is the angle-averaged value of the coefficient 
ci(d;ni), 

d (d) := <£ ci(d;ni)dfl d _i. (A28) 

Js d ~ 1 ( 0,1) 

One checks that always 

lim ci(d) = Ci(3) = ci. (A29) 

d—> 3 

The radial integral in Eq. (1A27I) is convergent if the real 
part 5ft(d) of d fulfills 5ft(d) < 3. Let us introduce 

e := d — 3 

and let us expand ci(d) around e = 0 , 

ci(d) = Ci(3 + s') = Ci + c' 1 (3)e -I- 0(£ 2 ). (A30) 


Then the expansion of the integral /i(d) around e = 0 
reads 

hid) = k(3) + Cl In ^ + 0(e). (A31) 

2s 2 Iq 

In Eqs. (IA30I) (IA31I) we have used (IA29D . Let us note 
that the coefficient c' : (3) usually depends on In ri 2 and it 
has the structure 

ci(3) = c' n (3) + ci 2 (3)ln^. (A32) 

to 

Therefore the DR correction also changes the terms oc 
In f 12 • 

The DR correction to the RH-regularized value of the 
integral / RH ( 3 ;ei,£ 2 ) relies on replacing this integral by 

/ rh ( 3; ei, e 2 ) + A/i + A/ 2 , (A33) 

where 

AJ 0 := 4(d)-7^(3; £1 ), a = 1 , 2 . (A34) 

Then one computes the double limit 

lim (/ RH (3; d, e 2 ) + A/i + AJ 2 ) = A - Cl + — 

ei—>-3 \ / ZS 

£-2 —^3 

~ 9 ( c i(3) + c' 2 (3)) + (ci + c 2 ) In + 0(e) 

Z to 

= ^- ~2e~ ~ 2 ^ + c ' 21 

+ ( C 1 - 2 C 12 ( 3 ) + c 2 - 2 C 22 ( 3 )) In + 0(e). (A35) 

Note that all poles oc 1/ei, l/e 2 and all terms depending 
on radii i\, li or scales Si, S 2 cancel each other. The 
result (IA35I) is as if all computations were fully done in 
d dimensions. 

In the DR correcting UV divergences of the 3PN two- 
point-mass Hamiltonian performed in Ref. [31] , after col¬ 
lecting all terms of the type (IA35I) together, all poles 
oc 1 /(d— 3) cancel each other. This is not the case for the 
UV divergences of the 4PN two-point-mass Hamiltonian 
considered in the present paper. As explained in Sec. 

IVIIIDI of our paper after collecting all terms of the type 

ASH) , one has to add to the Hamiltonian a unique total 
time derivative [given in Eq. (18.241) ] to eliminate all poles 
oc 1/ (d — 3) (together with 7 0 -dependent logarithms). 

3. IR corrections 

To regularize IR-divergent integrals which appear in 

the derivation of the 4PN two-point-mass Hamiltonian, 

we have originally developed a technique analogous to the 

one described above and used to compute DR corrections 

to UV-divergent integrals regularized in 3 dimensions. 

After completing tedious computations we have obtained 
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IR terms analogous to UV terms described by Eq. (1A35I) . 
After adding all these terms we have the expression with 
poles oc l/(d — 3). Then we have checked that there 
exists no total time derivative by means of which one can 
eliminate these poles. The conclusion was that even in 
d-dimensional computation of IR-divergent integrals one 
has to introduce an additional IR regularization factor 
( r/s) B with a new scale s. 

To be consistent with DR correction of UV diver¬ 
gences performed in d dimensions, we have developed 
a d-dimensional version of IR regularization. We have 
devised two different regularization schemes. We will 
however see that the results of these regularizations are 
identical with the results achieved by means of purely 
three-dimensional computations. 

There is a crucial difference between the results of ap¬ 
plication of UV and IR regularizations: the results of IR 
regularizations depend on an arbitrary IR regularization 
scale which we had to introduce, whereas the result of UV 
regularization does not depend on any scale. Moreover 
the results of two different IR regularizations developed 
below are different; therefore, we have to conclude that 
the result of IR regularization is ambiguous. In Appendix 
IA 3 cl we show that the ambiguity can be expressed in 
terms of only one unknown dimensionless parameter. 

The basic idea of both IR regularizations is, similarly 
to what we have done for UV divergences, to replace 
this part of the three-dimensional integral /^, H (3; a, b, e) 
from Eq. (lAdl) which is responsible for IR divergences, by 
its d-dimensional counterpart. Let us thus consider the 
three-dimensional integral (1A4I) which is IR divergent. 
In all considered cases one checks that its IR pole term 
proportional to 1 /e is related to this part of the integrand 
*(x) of dSS): which, after expansion around r = oo (r := 
|x|), is proportional to 1/r 3 , 


a. Modifying behavior of the d 3 at infinity 


All terms contributing to poles at spatial infinity are 
collected in Eq. (I7.12bl) . They (with the exception of the 
first term) have the structure 

(A39) 

and the first term in Eq. (17. 12bl) we treat as 
2 (d-1) ^) h wa^ A d 

To regularize all these terms properly we have made the 
replacement 


a- i7 TT 
ft (4 )ij 



(A40) 


where ( r/s) B is an IR regularization factor with s being 
a new constant (needed to make the regularization factor 
dimensionless). 

After making the replacement (1A40I) in each term in 
Eq. (I7.12bl) . we have found in d-dimensions for each 
term this part of its expansion around r = oo which 
contributes to the IR divergence. It is proportional to 
r 6-3 d+B an( j } las the structure 




-'TT 

\i )ij 


x s~ B r 6 ~ 3d+B 


+ ci(d,R;n) 
when r —> oo. (A41) 


One can make the replacement (IA40I) directly in d = 
3 dimensions, then instead of the expansion (1A41I) one 
obtains 


*(x) = • • • + Coo(n) r 3 -|-, when r —> oo, (A36) 

where n := x/r. It means that one can reproduce 
the pole term of Eq. (IA41) by means of integration of 
Coo(n)r~ 3 taken over the exterior of the ball B(0,i?) 
(with the center at the origin 0 of the coordinate sys¬ 
tem and the radius R). We thus define 

^S> H (3; a, b, e) := f Coo(n) r~ 3 (—V d 3 x 

Jr 3 \m(o,r) Vr o / 


—3 


= Coo r 

Jr 




‘dr 


= -c 0 


+ In —+ 0(e), 


(a + b)e r 0 ) 


(A37) 


where Coo is the angle-averaged value of the coefficient 

Coo(n), 


CqO • — 


j) Coo(n) dfU- (A38) 


s 2 (o,i) 


Aj(x) a 1 (J-^j 


■ B h TT 
ft (4 )ij 


x s~ B r B ~ 3 


+ c^(B; n) 
when r —> oo. (A42) 


We have checked that always 


lim c < ( 0 (d,R;n) = c^(R;n), (A43a) 

d—> 3 

lim c^(B; n) = Coo(n), (A43b) 

B—t 0 

where Coo(n) is the coefficient in the three-dimensional 
expansion (IA36I) . 

For all terms from Eq. (I7.12bl) we have computed the 
integral 

IL(d, B) := s~ B f c^(d,B-n) r 6 ~ 3d+B d d x 

jR d \B(O.R) 


r c 

= c 1 oo(d,B)s- B 

Jr 


.5-2 d+B 


dr 


= -<4, {d,B)s 


-B 


^6-2 d+B 

6-2 d + B' 


(A44) 
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where c^ 0 (d, B) is the angle-averaged value of the coeffi¬ 
cient cl(d,B;n), 

4 (d,B):=(f c^(d,B;n)dn d _i. (A45) 
Js d ~ 1 (o,i) 

One easily checks that the equalities analogous to (1A43I) 
are also fulfilled for the angle-averaged values of the co¬ 
efficients, 

]im 4 (d, B) = 4 (B) , (A46a) 

d—>3 

Um 4 (B) = Coo , (A46b) 

where the three-dimensional coefficient 4 (-®) is the 
angle-averaged value of the three-dimensional coefficient 
c£CB;n), 

44):=/ 4(B;n)dfi 2 . (A47) 

ds 2 ( 0 ,l) 

The integral 4(d, B) can be shortly written as (let us 
recall that e •= d— 3) 

I 1 00 (d,B) = I 1 00 (Z + e,B) = (A48) 


In view of the equalities (IA46I) it is clear that the value 
of the right-hand side of Eq. (IA52I) would be the same 
if all computations leading to it were performed in 3 di¬ 
mensions. We thus have 

[ i(x) d 3 x^j = A-^^-( 0 )-c oo ln—. (A53) 
Jr* J reg 1 dB s 

b. d-dimensional RH regularization 

We have also considered another way of regularizng IR 
divergences. Namely, before integrating an IR-divergent 
integral over d-dimensional space, we multiply the full 
integrand by a factor 

(?n?r <-> 

and after evaluating it we take the finite part of the IR 
pole occurring at a\ + a-i = 2(d — 3). This recipe means 
nothing more than the d-dimensional version of the Riesz- 
implemented Hadamard regularization we performed in 
3 dimensions. 

In this approach instead of the expansion (IA41I) one 
has (we introduce here /3 := aq + a^) 


where we have defined 

N(e, B) ■= -<4,(3 + e, B)s~ B R B - 2e . (A49) 

As the regularized value of the integral (1A48I) for e —> 
0 and B —> 0 we take the finite part (FP) of the pole 
occurring at B = 2e in d dimensions (we follow here Ref. 
[60i] : see especially Sec. VIII there). We thus define 


FP/i, := FP linr 

e-s-0 B-yO 


N(e,B)-N(e, 2e) 
B- 2e 


Fr N(e, 0) — N(e, 2s) 
e —>-0 —2 s 


= FP 

£—>-0 


dN 


—( 0 , 0 )+ 0(4 =—( 0 , 0 ) 


dB 


dN, 


dB 


/? 

-gfM ~4 (3,0) In-. 


(A50) 


/ij( x ) A d = '' • + 4(4 n) s 0 r 6 3d+/3 + • • • , 

when r —► oo, (A55) 

and instead of the integral (IA44I) one considers the inte¬ 
gral 

4 (d, 0) := s-P f 4 (d; n) r 6 ~ 3d +$ d d x 

jR d \M(0,R ) 

/»oo 

= 4 (d)s"M r 5 ~ 2d +P dr 

Jr 

o6-2d+/3 

= -^*-VanHr (A56) 

where 

4(4 := / 4( d ;n)dfi d -i- (A57) 

ds^-Rop) 


As the correction to integral / RH (3; a, 6 , e) we define the 
difference 


A 4 := FP 4 — / RH (3; a, b, e), (A51) 

so the regularized value of the three-dimensional IR- 
divergent integral over *(x) reads 


[ i(x) d 3 ^ := lim (d RH (3; a, b, e ) + All) 

V y regl 


= A- -gfM -4(3,0)In -f. 

(A52) 


One checks that always 


lim c^(d) = 4,(3) = Coo 


(A58) 


The crucial difference between the integral 4(d, 8) 
and the integral ll(d,B) of (IA44D is such that in 
4 (d, 3) the coefficient 4(4 n ) [and its angle-averaged 
value 4(d)] does not depend on /?, whereas in Il(d,B) 
the coefficient 4 (d,B;n) [and its angle-averaged value 
4(d;.B)] does depend on B. The integral (1A56I) can be 
written as 


4 ( 4 / 3 ) 


p-2e’ 


(A59a) 
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ri(e,/3) := -c^(3 + e)s 2e , (A59b) 


and its regularized value for e —> 0 and /3 —» 0 we define 
as the finite part 


FP / 2 := FP lim 

00 e-X) /3-X> 


?7(e,/3) - rj(e, 2s) 
P-2e 


= l?( 0 ’ 0 ) = - c -^ 3)ln -- ( A60 ) 

op s 

The correction to the integral / RH (3;a, 6, e) is again the 
difference 

All ■■= FP II - /™(3; a, b 7 e), (A61) 

so the regularized value of the three-dimensional IR- 
divergent integral over i(x) reads 

(f ifx.)d 3 x] := lim (/ RH (3; a, 6, e) + All) 

\J R 3 / reg 2 £ ^° 

= A- 4(3) In r -f. (A62) 

By virtue of Eq. (IA58I) it is clear that the value of the 
right-hand side of Eq. (1A63I) would be the same if all 
computations were performed in 3 dimensions. We thus 
have 

f f i(x)d 3 x] = A — Cqo In(A63) 

\J R 3 / reg 2 S 

Let us finally note that the above result can be immedi¬ 
ately read off from Eq. (IA4I) after dropping the pole part 
and identifying the scales r$ and s. 


c. IR ambiguity 


4. “Partie flnie” value of singular function 


The concept of “partie finie” value of function at its 
singular point was previously used e.g. in the calcula¬ 
tion of ADM point-particle Hamiltonians at the 2PN and 
2.5PN levels (see Appendix B in 75}), and also at the 
3PN and 3.5PN orders (see also (ilfh 

Let / be a smooth real-valued function defined in an 
open ball B(x 0 ,E) C R 3 of radius E > 0 and origin at 
x 0 € R 3 , excluding the point x 0 , i.e., / € C' 00 (B(xo, E) \ 
{x 0 }). At x 0 the function / is assumed to be singular. It 
is enough to consider functions / which have only rational 
singularities. We thus assume that there exists positive 
integer N such that the limit 

lim /(x)|x — xol^ (A 66 ) 

x->x» 

exists and is finite. Let us denote by N m - m the smallest 
positive integer N for which the limit (IA 66 I) exists and is 
finite. We define 


fl(x) 


f (x) |x X 0 | Mnin ^ for X^X 0 , 

lim x _> Xo / (x) | x - x 0 1 ■ Nmin , for x = x 0 . 


(A67) 

We also define two families of auxiliary functions / n and 
g n (labeled by unit vectors n, n • n = 1 ): 


(0; E) 9 e / n (e) := / (x 0 + en) <E R, (A 68 a) 

(0; E)B e g n (e) := g (x 0 + en) e R. (A 68 b) 


Our final assumption is that for all unit vectors n the 
function g n is smooth in the half-closed interval (0;E). 
Using Eqs. (IA67D and (IA 68 I) we obtain 


Comparison of the two IR regularization schemes con¬ 
sidered above immediately leads to the conclusion that 
the results of their application to computation of the in¬ 
tegral of h^\ 2j ■ Eq. (17.12bl) . are different. By virtue of 
Eqs. (IA53I) and (IA63I) one gets 


Aff r eg2,2 

AA - n 4PN 



= -E|f(o). 


(A64) 


where the summation is over all terms of the integrand 
h^' l2 y We have computed the difference . It 

turns out that it can be written as 


137 

Af/^p 8 ^’ 2 = (total time derivative) + ——F, (A65) 


768 


where F is defined in Eq. (18.141) . 

To take into account ambiguity (1A65I) of IR regular¬ 
ization we have introduced in Eq. (18.111) a dimensionless 
ambiguity parameter C. According to Eq. (IA65I) this am¬ 
biguity can be written (up to adding a total time deriva¬ 
tive) as a multiple of the term F, and this is the content 
of Eq. (18.181) . 


/n(e) = (A69) 

£4 v min 

From Eq. (1A69I) . after expanding g n into Taylor series 
around e = 0, the formal expansion of / n into Laurent 
series around e = 0 follows: 


/n(e) = J2 


(n), 




ffn(O) ffn(O) , , 9 

AT I AT 1 I ’ ’ ’ —I 


(Afmin-l) 


( 0 ) 


(7V min -l)!, 


1 


N ■ > 

1 y min • 


gi N ™»\0) + O{e). 


(A70) 


We define the regularized “partie finie” value of the func¬ 
tion / at Xo as the coefficient at e° in the expansion (1A70I) 
averaged over all unit vectors n: 

/reg(x 0 ) := ^ dU 2 a 0 (n) 

= 3ski/ dn29 ”""‘" >(0) ' 


(A71) 
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Let us consider the function 

g(x) := f(x)\x-x 0 \ N , 

where N is the positive integer N > N m - ln Then one can 
define / reg (xo) using the function g instead of g. It is 
easy to show that the value of / reg (xo) will not change, 
so the definition (IA71I) does not depend on the choice of 
the number N provided N > -/V m ; n . 

Definition (1A71I) is used to give the meaning to inte¬ 
grals of /(x) S a . where the function / is assumed to be 
singular at x = x a . Namely, we define 

J d 3 x /(x) Sa ■■= /reg(x a ). (A72) 

The definition (IA72I) is an extension of the notion 
of “good” Dirac 5-functions introduced by Infeld and 
Plebanski [jdj . Their definition assumes that “good” S- 
function, besides having the properties of ordinary Dirac 
S distributions, also satisfies the condition (cf. Appendix 
1 in p4j ) 

d 3 a; = 0, for k = 1, 2 ,.. . ,p. (A73) 

^ a 

Obviously the definition (IA72I) entails the fulfillment of 
the condition (IA73D . 

The important feature of the definition (IA71I) is that 
the regularized value of the product of functions is not, 
in general, equal to the product of the regularized values 
of the individual functions: 

(/ 1/2 • • • ) reg (x 0 ) ^ fi reg (xo) h reg (xo)-" (A74) 

(the above property but with the equality sign was called 
“tweedling of products” by Infeld and Plebanski; see, 
e.g., Appendix 2 in 0 ). The property (IA74I) of three- 
dimensional “partie finie” operation leads to problems. 
Let us consider some function S singular at x = x a . It 
is natural to demand that 

S(x) S a = 5 reg (x 0 ) 5 a . (A75) 

The above rule is always used when solving Poisson equa¬ 
tions with singular source terms of the form S (x) S a [see 
Eq. (1B2D below]. Then, multiplying both sides of Eq. 
(IA75I) by another function T which is also singular at 
x = x a , one gets 

T(x) S(x) S a = T(x) S' reg (x a ) S a - (A76) 
The rule (1A75D applied to Eq. (IA76D implies that 


generic d-dimensional case the distributivity is always 
satisfied, 

(/l d) /i d) • • • ) reg ( X a) = /l ( reg( X a)/2 ( reg(Xa) ‘ , (A78) 

where /reg (x a ) is rather not defined by the d-dimensional 
analog of the definition (IA711> , but it is directly computed 
by employing the existence of such region in the complex 
d-plane where the function /^ (x) is finite for x = x Q (see 
the example at the end of this appendix). In d dimensions 
one can thus always use 

/ (d) (x) 6^ d \x - x„) = /W(x a ) S^ d \x - x a ). (A79) 

Therefore one defines the dimensional-regularization 
rule, 

/reg(x a ) := lim (/^ (x a )) , (A80) 

where /^ is the d-dimensional version of /. In the com¬ 
putation of the 4PN Hamiltonian the usage of the rule 
(IA80I) boils down in practice to the usage of the three- 
dimensional definition (1A71I) after the usage of distribu¬ 
tivity (IA78II (in the case of computing regularized value 
of the product of functions). 

As an example of justifying the distributivity (IA78I) 
let us consider the following 4PN-related contact integral 
[here we denote by < f>^] the function, given by the right- 
hand side of Eq. (ICla[l , which is the d-dimensional ver¬ 
sion of the three-dimensional potential denoted by 4 >( 2 )]'- 

J A d x{^j>^ (x)^ 5i = J d d x K 5 (mir'l~ d + m 2 rl~ d )° Si. 

(A81) 

Because 5f(d) < 2 => lim x _> Xl d = 0, then 

J d d x(^ 2 )( x )) 5i =« 5 (m 2 r^ d ) 5 = [(^^^(xi)] 5 . 

(A82) 

Therefore the three-dimensional value of the integral 
(IA81[) equals 

J d 3 x (^( 2 ) (x )) 5 5, = lim [(^ 2 j) reg (xi)]° 

= [(</>( 2 )) reg (xi)]°- (A83) 

Computation of the three-dimensional integral of 
( 0 ( 2 ) (x)) 5 5i directly by means of Eq. (1A71I) leads to the 
result different from the result of Eq. (1A83I) . 


(TS) leg (x a ) =T leg (x a )S leg (x a ), (A77) 5 - Distributional differentiation of homogeneous 

functions 

what in general contradicts Eq. (IA74I) . 

To avoid this kind of problem one should employ di- Appearance of UV divergences is not the only cost of 
mensional regularization. One can check that in the employing Dirac-delta sources. Another consequence is 
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that in our computations we have to differentiate homo¬ 
geneous functions using an enhanced (or distributional) 
rule, which comes from standard distribution theory (see 
Sec. 3.3 in Chapter III of Ref. @)- 

Let / be a real-valued function defined in a neighbor¬ 
hood of the origin of R 3 . / is said to be a positively ho¬ 
mogeneous function of degree A, if for any number a > 0 

/(ax) = a A /(x). (A84) 

Let k := —A — 2. If A is an integer and if A < —2 (i.e. k 
is a non-negative integer), then the partial derivative of 
/ with respect to the coordinate x l has to be calculated 
by means of the formula [cf. Eq. (5.15) in 0 ] 


<%/(x) 


<V(x) + 


(-l)fc d k S( x) 

k\ dx 11 ■ ■ ■ dx lk 


As an example let us employ the formula (IA89I) to cal¬ 
culate first and second partial derivatives of l/r a , l/r^, 
and 1/r 3 . For the first partial derivatives we obtain 


1 

di-, 

(A90a) 

r a 

~r a 

1 

fa — 1 

' a 

T a 

(A90b) 

1 

1 4:7 r _ 

(A90c) 

Z 3 " 1 

' a 

di 77 - —d t d a . 

~ r l 3 


Let us note that in Eq. (lA90a,H there is no need to use 
the rule (IA89I) . and in Eq. (IA90bl) the term proportional 
to d a [obtained from the usage of (IA89I) ]. vanishes. The 
second partial derivatives read 


x d m /(x') x' n ■ ■ ■ x ,lk , (A85) 


did,- = d^- - ^M«, 

r a ~ J -r a 3 


(A91a) 


where dif means the derivative computed using the stan¬ 
dard rules of differentiations, E is any smooth close sur¬ 
face surrounding the origin and deq is the surface element 
on E. 

The rule (IA85I) should also be applied to differentiation 
of functions which are homogeneous not with respect to 
x, but with respect to x-xo, for some constant Xo G R 3 . 
Let / be such function, then there exists another function 
cA for which the relation 

/(x) = 0(x - x 0 ) (A 86 ) 

is fulfilled and the function £ >->■ (/)(£) is a positively homo¬ 
geneous function of degree A, i.e., for any number a > 0 

</>(£)■ ( A87 ) 


d < d if2 = ( A91b ) 

'a 'a 

didj = didj -4 - ^ (1 6didjS a + 3dijAd a ). (A91c) 

r a r a 

Making use of Eqs. (IA91al) . (IA91bl) . and (lA91cl) . one gets 

A— = — 47 t <5 a , (A92a) 

r a 

^ (A92b) 

'a 'a 

.1 6 107T 

A- = - - — A8 a . (A92c) 


Obviously (here £ := x — x 0 ) 

df{yj) = d(f(£) 
dx 1 df l 


(A 88 ) 


and for the derivative dcf/df 1 the rule (IA84I) is directly 
applicable. Therefore the partial derivative with respect 
to the coordinate x l of the function / satisfying condi¬ 
tions (1A86I) and (1A87I) should, by virtue of (1A88I) and 
(IA85I) . be calculated by means of the formula 


diffx) = dif (x) + 


(—l) fc d k d(x — x 0 ) 
k\ dx 11 ■ • ■ dx lk 


x£d a^(Oe ,ll ---r fc , (A89) 


where d^f means the derivative computed using the stan¬ 
dard rules of differentiations, E is any smooth close sur¬ 
face surrounding the point Xo and dcq is the surface ele¬ 
ment on E. 


The distributional derivative (1A89II does not obey the 
Leibniz rule. It can easily be seen by considering the 
distributional partial derivative of the product 1 /r a and 
l/r„. Let us suppose that the Leibniz rule is applicable 
here: 


di -7T - di 



1 ~ 1 ^ 1 , 1 

—y di -1- di — . 

r a r a r a ri 


(A93) 


By virtue of Eqs. (IA90al) and (IA90bl) the right-hand side 
of Eq. (1A93I) can be computed using standard differential 
calculus (no terms with Dirac deltas), whereas comput¬ 
ing the left-hand side of (1A93I) by means of (lA90cl) one 
obtains some term proportional to did a - 

The distributional differentiation is necessary when 
one differentiates homogeneous functions under the inte¬ 
gral sign. Let us consider the following locally divergent 
integral (here a yf b): 


Pai Paj 


d 3 x 



(A94) 
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We shall regularize this integral in two different ways. 
We first replace in (IA94I) differentiations with respect to 
x l by those with respect to x l a (obviously diT a = — d a ir a )- 
Then we shift the differentiations before the integral sign 
and apply directly the Riesz formula (IA5I) . The result is 


PaiPaj I d xfdidj 'j 4 — PaiPaj daidaj [ 4 

J \ ' o, J 'b J a b 


Pai Paj C'aiC'aj 


i daifiai I 


27T 


ab 


4?r[p 2 - 4(n ah • p Q ,) 2 ] 


(A95) 


' ab 


We have obtained (1A95I) performing integration first and 
then differentiation. Now we shall regularize the integral 
(IA94I) doing differentiation first. To do it we have to use 
the rule (IA89I) . which gives [cf. Eq. (lA91al) ] 

didj ^ = (3 nini - S tj ) 1 ^ S ZJ S a . (A96) 

I a > n 'A 


We substitute (IA96I) into (IA94I) : 


Pai Paj 


d 3 :r 



1 


r 


4 

b 



3(n Q • p a ) 2 - p 2 
q 4 



d 3 x^4. (A97) 

r b 


^d -2 — ®i-®L r d-2 nd^ i ^ a ' (AlOOb) 

Using the formula (1A85I) [or (1A89D ] in d (or 3) di¬ 
mensions requires averaging of products of unit vectors 
over the unit sphere. This can be done by means of 
the following formulas (see Appendix A 2 of 28] for d- 
dimensional expression and Appendix A 5 of [80| for its 
three-dimensional version): 


dU d _i n h ■ ■ ■ n i2k+1 =0, fc = 0,1,2 ,... , 

(AlOla) 

-i nil ' ‘ ‘ nl2k = (d + 2(k- l))!!^” 1 

x ^{* 1*2 ‘'' ^i2k-ii2k }> A: — 0, lj 2,... . 

(AlOlb) 


Here, for positive integer k, fc!! := k(k — 2) • • • 1 for k odd 
and fc!! := k(k — 2) • • • 2 for k even, Sl^-i is the area of 
the unit sphere in 


l = 


27T d / 2 

W’ 


(A102) 


and A {ili2 ... ik } := J2aes A <r(ii)-<r(i h ), where S is 
the smallest set of permutations of {l,...,fc} making 
A{hi 2 ---ik} fully symmetric in Ufc- Let us note 

some simple cases: 


The second integral on the right-hand side of (1A97I) is 
calculated by means of the definition (1A72I) . The obvious 
result is 

= ( A98 ) 

j 'b 'ab 

To calculate the first integral on the right-hand side of 
(1A97I) we apply the procedure described in Appendix I A II 
We obtain 

f j3 3(n a • p a ) 2 - p 2 16 tt [p 2 - 3(n ah • p a ) 2 ] 

a X T .3 r 4 _ o„4 

4 'a'b °'ab 

(A99) 

Collecting Eqs. (IA97I) (IA99I) together we get the result 
which coincides with the result (1A95I) obtained before. 
The two ways of regularizing the integral (1A94I) . de¬ 
scribed above, coincide only if we apply formula (IA89I) 
when we perform differentiation before integration. 

It is not difficult to show that the formula (1A85I) is also 
valid (without any change) in the d-dimensional case, i.e. 
when / is a real-valued and positively homogeneous func¬ 
tion of degree A defined in a neighborhood of the origin 
of R d . The formula (1A85I) can be applied in d dimensions 
when the number fc:=—A + d—lisa non-negative in¬ 
teger. For example, the d-dimensional versions of Eqs. 
(lA90al) and (lA91al) read 

di-jz 2 = ’ (AlOOa) 

T a Ta 


d 


(A103a) 


dfld-in ll n ,2 n l3 n li = 




d -1 


d(d + 2) 

^ (^1^2^23*4 T dq d~ Si 1 i 4 di 2 i 3 '). 

(Al03b) 


6. Riesz kernel 


To avoid problems related with applying distributional 
derivatives (even in d dimensions) one can replace Dirac 
delta distribution S by some functional representation 
(“delta sequence”). In d dimensions one can e.g. employ 
the Riesz kernel K a (e a )- 


S a = lim K a (e a ), K a (e a ) 
£ a —>• 0 


r((d- e „)/2) £a 

7r d / 2 2^T(e a /2) ra 


(A104) 

Then one should replace in the constraint equations m 
Dirac-delta sources by Riesz kernels (1A104I) . solve the 
constraints perturbatively and develop the whole PN 
scheme (let us stress that then no distributional differ¬ 
entiation is needed). At the end of the calculation, one 
takes the limits ei —> 0 , 62 —> 0 , and only after this one 
computes d —> 3 limit. 

This procedure was applied by us in [4(| to recompute 
all UV divergent integrals at the 3PN level. It is however 
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too difficult to be performed fully at the 4PN level, so it 
has not been applied in the main part of the present pa¬ 
per, but it was used for checking some results. It should 
be mentioned that the Riesz-kernel method of regular¬ 
ization has been applied by Damour in Refs. 1781 [79j. It 
also may be pointed out that_ the dimensional regular¬ 
ization calculations of Ref. [63] have been performed in 
momentum representation, where also only ordinary (i.e., 
nondistributional) space-time derivatives show up. 

As an example let us first compute the potential (f)( 2 ) 
for Riesz-kernel sources. Instead of Eq. (16. 10 all we thus 
solve 

A d 0( 2 ) = -^rriaKa. (A105) 

a 

Making use of Eq. (1B41) one gets 


Appendix B: Inverse Laplacians 

In the present paper we have to consider (both in d 
and in 3 dimensions) numerous Poisson equations with 
distributional source terms, 

A df = E«( X )^ ( B1 ) 

a 

where usually the function g is singular at x = x a . Equa¬ 
tion of this type we solve as follows: 

/( X ) = Ad 1 5a = A d 1 J29reg(Xa) S a 

= Sreg(x a )A rf S a , (B2) 


h 2 ) 


r r “ *- d/2 m-ea)/2) 2 _ d+e , 

V £a (2 — d + e a )r(e a /2) a “ 


(A106) 

Let us next consider the integral of Eq. (IA81I) , which now 
takes the form 


Jd^x^Ki. (A107) 


where the regularized value <? reg of the function g is de¬ 
fined in Eq. (IA71I) (in 3 dimensions) or in Eq. (IA80I) 
(when, to avoid ambiguities, one has to employ the d- 
dimensional version of computing of g le g)- The function 
A^ 1 ^ is defined below in Eq. (IB9I) (in d dimensions) or 
in Eq. (fBTSl) (in 3 dimensions). 


To compute this integral it is enough to use Eq. dMt - 
After computing the double limit e± —> 0, 62 —> 0 (the 
order of the limits does not matter) one obtains 


1. d-dimensional inverse Laplacians 

We start from considering some solutions of equation 


7r- sd / 2 r(d/2) 5 

32(d — 2 ) 5 


5„10-5d 


m 2 r 12 


(A108) 


The value of the above formula in the limit d —> 3 coin¬ 
cides with the result (1A83I) . 

Let us finally mention that the usage of the Riesz kernel 
directly in 3 dimensions does not resolve ambiguities. For 
example, three-dimensional computations corresponding 
to d-dimensional ones presented in the above example 
would lead to the result which is different from (1A83I) 
(so the “tweedling property” of the product would be 
violated). 

I 


A 3/ = r“, n= 1 , 2 ,..., (B3) 

where A d denotes n-folded superposition of the d- 
dimensional flat Laplacian A^. By direct computation 
one checks that the function 

A~n a , = r(a /2 + l)r((a + d)/ 2 ) 2n 

d a ' 2 2 "r(a/2 + n + l)r((a + d)/2 + n) a 

(B4) 

is the solution of Eq. m- The function A d n r% we call 
nth inverse Laplacian of r“. Let us note the formulas for 
the first, second, and third inverse Laplacians of r“: 


A 

A 

A 


-v . = 

d a ’ 


- 2 r a 
d 'a ■ — 


-3 r a 
d a 


_ _a _ 

(of + 2)(a + d) 

r 4+o: 

_ _a _ 

(2 -\- a)(4 -|- ql)(cL + a)(2 d o;) 

r 6+a 

_ _a _ 

(2 a)(A q )(6 + OL){d H- gi)(2 H- d c^)(4 d -\~ c^) 


We also have to consider solutions of equation 


A 3/ = s t 


n= 1 , 2 ,... . 


Making use of Eq. (IB41) and the formula 


(B5a) 

(B5b) 

(B5c) 


(B 6 ) 


A d rl d = -k 1 d, 


(B7) 
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one checks that the function 


A7"«5 q := - 


nY{2-d/2) 


2n—d 


4 n_1 (n — 1)! T(n + 1 — d/2) a 


(B 8 ) 


We also need solutions of equations of a more general 
type than Eq. (ESI). It reads 

A 3/= r>£ ■■■«?, n= 1,2,.... (BIO) 


solves Eq. (IB 6 I) . The function At"< 5 a we call the nth 
inverse Laplacian of S a - The first inverse Laplacian of S a 
thus reads 

A ^S a :=-nr 2 a - d . (B9) 

I 


More precisely, we need to consider Eq. mm for n = 
1,2,3 and for each of these values of n we need to take k = 
1,... , 6 . To save space we list below solutions (inverse 
Laplacians) related to Eq. (IB10I) only for n = 1, 2,3 and 
k = 1,2, 3: 


r 3 r a 

L a' a 


r 2+ot 

_ "a 1 a _ 

(ex + l)(oi + d + 1) 

[(2 + a)(d + a)n l a nj - 2 %]rg+ a 
ol(2 -|- q)(c? + a)(2 -\- d - 1- O') 


<~ 1 n i n j n k r a 
,l a n a n a'a 


[ 2 (^fc< + 5 lk n{ + 8 jjn k a ) - (1 + a)(l + d + o»>q]rg +a 
(1 — cy)(l “I - c^)(l d a)(3 -\- d -\- oi) 




A„, 2 n!, nii 


V ^77* 71 ^ 71 ^ V°^ 


r 4 +a 

__ "a a _ 

(1 + cy)(3 + c^)(l H - d a)(3 H- d cy) 

[(4 + a)(d + a)n> 3 a - 4^-]r^ +a 
a (2 -|- q)(4 + ol){(L + a) (2 d cy)(4 -|- d cy) 

[(3 + a)(l + d + a)n»k - 4(^-n^ + <5j fc n» + ^fcn^)]r^ +a 
(q — 1)(1 T c^) (3 T q)(l T d T q)(3 T d T q)(5 T d T q) 




A, 3 n l „n 3 „i 


^d 


r 6+Oi 

_ "a 1 a _ 

(1 + cy)(3 + q)( 5 + cy)(l ~\~ d -\- a)(3 d q)(5 d cn) 

_[(6 + a)(d + ct)n>4 - 6^]rg+ a _ 

a(2 H - a)(4 -|- cy)(6 H- cn)((i -|- cy)(2 -|- d H- g:)(4 -|- d -b q)(6 d -\- o;) 

_[(5 + q)(l 4- d + a)n l a n 3 a n k a - 6(i^n* + S jk n l a + &fcn£)]rg+ a _ 

(q — 1)(1 4~ rr)(3 4“ q)(5 T q)(l T d T q)(3 4- d 4- q)(5 -t- d T q)(7 T d T q) 


(Blla) 

(Bllb) 

(Bile) 

(Blld) 

(Bile) 

(Bllf) 

(Bllg) 

(Bllh) 

(Bill) 


2. Three-dimensional inverse Laplacians 


The function 


For convenience let us start from rewriting some d- 
dimensional results of Appendix IB II in d = 3 dimensions. 
The function 


A -ra a ._ T(q + 2) g+2n 

a ‘ T(q + 2 n 4 - 2 ) a 


(B12) 


A _n (5 a := - 


2n—3 


47 r( 2 n- 2 )l a 


(B14) 


solves Eq. flBfill in d = 3 dimensions and is the nth inverse 
Laplacian of 5 a ■ The first inverse Laplacian of 5 a thus 
equals 


is the nth inverse Laplacian of r“ [it is thus the solution 
of Eq. (IB3I) in d = 3 dimensions], so the first inverse 
Laplacian of r“ reads 


rpOt~\~ 2 

A _ 1 r“ -- — - 

a (a + 2) (a + 3) 


(B13) 


A - 1 (5 a := -- 2 —. (B15) 

47rr a 

In the derivations of the field functions needed to cal¬ 
culate the 4PN-accurate two-point-mass Hamiltonian we 
have used some special solutions to the partial differential 
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equations of the form 

A n f = r a y b (a^h), n = 1 , 2 ,.... (B16) 

This special solution to Eq. (IB16I) is defined below, we 
denote it by 7(n; a, (3) and call the nth inverse Laplacian 

°f r a r b ' 

Let us mention that in Ref. [81| one can find formulas 
[see Eq. (B13) there and equations below it], which de¬ 
scribe the formal solution to Eq. flBTH) in d dimensions 
for any complex values of n, a, and /3. The solution is 
expressed in terms of the Appell hypergeometric func¬ 
tion F 4 of two variables. We have tried to use it in our 
4PN-related computations but, in general, this has led to 
calculations too complicated to be useful. 

The large family of three-dimensional inverse Lapla- 
cians of r“r^ is based on the solution to the Poisson 
equation 

A / = —. (B17) 

r a n 

One immediately checks that ln(r 0 + r b ± r a b ) solves Eq. 
m LB- Because ln(r a + rb — r a b) is singular along the 
segment joining the points x Q and x b , the first inverse 
Laplacian of r“ 1 rj[~ 1 we define to be 

7(1; -1, -1) := ln(r a + r b + r ab ). (B18) 

Let us apply the operator A a (A a contains differenti¬ 
ations with respect to x l a only) to the both sides of Eq. 
mm- Assuming commutativity of the operators A a and 
A one obtains 

A a (A/) = A(A a /) = — (A a —^ = -47T—, (B19) 

r b \ r a ) r b 

where we have used 

A a — = A— = -4TrS a . (B20) 

r a r a 

After applying the rule [see Eq. (IA75I) and the discussion 
around this equation] 

f(x)5 a = /reg(x Q )5 a (B21) 

to the right-hand side of Eq. (iBTflll . one gets 

A(A a f) = -4t r^A. (B22) 

Tab 

Using Eq. (IB 1511 from (IB22I) one obtains 

A a f = -—A~ 1 6 a = —. (B23) 

Tab r a r ab 

Analogously one can derive the equation 

A&/ = —-—. (B24) 

r b r a b 


One checks that the solution (IB18I) fulfills Eqs. (IB23I) and 
(IB241) . whereas the function ln(r a + r b — r a b ) does not. 

Using the result (lBl8l) it is possible to calculate all 
multiple inverse Laplacians of the typ^B 

7(n; 2k — 1,21 — 1), n = 1,2,, k,l = 0,1, 2,... . 

(B25) 

To obtain I(n; 2k —1,2 1 — 1) (for n = 1,2,..., k, l = 
0,1, 2,...) one assumes that it has the following struc¬ 
ture: 

7(n; 2k -1,2 1- 1) := Wi(r a ,r b ,r ab ) 

+ W 2 (r a , r b , r ab ) ln(r a + r b + r ab ), 

(B26) 

where W\ and W 2 are polynomials of variables r a , rb, 
and r a b, consisting of terms of only 2 (k + l + n— l)-order 
in these variables. The form of Eq. (IB26I) can be inferred 
from dimensional analysis which takes into account the 
first inverse Laplacian of r~ 1 r given in Eq. (IB 181) . To 
fix the coefficients of W\ and W 2 uniquely one requires 
that the following conditions are fulfilled: 

AJ(n; 2k -1,21-1)= 7(n - 1; 2 k- 1, 21 - 1), 

(B27a) 

A a I(n-, 2k -1,21-1) = 2k(2k - l)/(n; 2k - 3, 22 - 1), 

(B27b) 

A b I(n-, 2k -1,21-1) = 21(21 - 1 )J(n; 2k -1,2 1- 3). 

(B27c) 

Equation (lB27bl) and (lB27cl) corresponds to Eq. (IB23I) 
and (IB24I) . respectively. In the case of k = 0 Eq. (IB27bl) 
should be replaced by 

r 2l ~ l 

A a I(n--1,21-1)= (27j ab _ 2) , ^ n - 3 ; (B27d) 

analogously for l = 0 instead of Eq. (IB27cl) one uses 

r 2fc-i 

A 6 7(n;2fc-1,—1)= a6 _ 2) , r^- 3 . (B27e) 

To illustrate the above procedure let us consider the 
second inverse Laplacian 7(2; —1, —1). We are thus look¬ 
ing for a solution / of the partial differential equation 
A 2 / = l/(r Q rb). We assume that the solution / is of 
the form W 1 (r a ,r b , r ab ) + W 2 (r a ,r b , r ab ) ln(r a + r b + r ob ), 
where W\ and W 2 are polynomials of indicated vari¬ 
ables consisting of only quadratic terms in these vari¬ 
ables (i.e. the most general polynomial of this type is 
air 2 + a 2 rl + a 3 r 2 h + a 4 r a r 6 + a 5 r a r ab + a 6 r 6 r ah ). We 


4 The method presented below was devised in Ref. H; see Ap¬ 
pendix C there. 
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need to fix the values of 12 coefficients (in fact less be¬ 
cause of symmetry with respect to interchanging the la¬ 
bels a and b of the particles) defining the polynomials W\ 
and W 2 . To do this we employ Eqs. (lB27al) . (IB27dl) . and 
(IB27cl) . which take the form: (1) A/ = ln(r Q + r b + r ob ); 
(2) A a f = r 0 /(2r ab ); (3) A b f = r b /(2r ab ). The equa¬ 


tions 1-3 fix the 12 coefficients of the polynomials W± 
and W -2 uniquely. The result is given in Eq. (IB28bf) be¬ 
low. 

Below we list explicit formulas for the inverse Lapla- 
cians which have been used throughout this paper (here 
a^b and s ab ■= r a + r b + r ob ): 


,-i 


1 


,-2 


l —3 


r a r b 

1 

r a r b 

1 

r a r b 


■— In S a bi 


(B28a) 


(~r 2 a + 3 r a r ab + r 2 ab - 3 r a r b + 3 r ab r b 
oo 


r b) + T^( r a~ r l b +r b )lns ab , 


(B28b) 


( - 63r„ + 150r 3 r ab + 126r 2 r 2 b - 90r a r 3 b - 63 r A ab - 90r 3 r b + 90r 2 r ab r b + 90r Q r 2 b r b - 90r 3 b r b 
126r 2 b r b - 90r a r b + 150r ab r b - 63r b ) 


A^ 1 — 

r b 

A" 2 — 
r b 


- 2r 2 r b + 90r a r ab r b 

+ ^ (3^ - 6r 2 r 2 b + 3r„ b + 2 r 2 a r 2 b - 6 r 2 ab rl + 3r£) In s ab , 

= {~rl - 3 r a r ab - r 2 ab + 3 r a r b + 3 r ab r b + r%) + ^ (r 2 + r 2 ab - r 2 b ) Ins ab , 

= (~ri - 30 r 3 a r ab - 62r 2 r 2 b + 90r a r 3 b + 63r 4 ab - 30r 3 r b + 30r 2 r ab r b 

- 90r a r 2 b r b + 90r 3 b r b - 62r 2 r b - 90r a r ab r b - 126r 2 b r b + 90r a r b 
+90r ob r b + 63r b ) + ^ ( r a + 2r 2 r 2 b - 3 r 4 ab + 2r 2 r£ + 6 r 2 ab rl - 3 r$) lns ab , 


(B28c) 

(B28d) 


(B28e) 


A ~ 3 — 
r b 


28^400 ( - 37r a - 210r 3 r a & - 951r^ 2 b + 1540r 3 r 3 b + 2013r 2 ^ b - 1050r a r 3 b - 1025r 3 b - 210r 3 r b 

- 210rfr ab r b - 840r 3 r 2 b r b + 840r 2 r 3 b r b + 1050r a r^ b r b - 1050r 3 b r b - 37'r*rg - 420r 3 r ab rf - 1062r 2 r 2 b rj 

- 2100r a r 3 b r b + 3075r 4 ab r 2 - 280r 3 r 3 + 420r 2 r ab r 3 - 2100r Q r 2 b r 3 + 2800r 3 b r 3 - 951r 2 r^ - 1050r a r ab r£ 

- 3075r 2 b r b + 1050r a rf + 1050r ab rf + 1025rf) + + 3r^r 2 b - 9r 2 r^ b + 5r® b + r 4 a rl + 6 r 2 a r 2 ab r 2 b 

~ 13r^ b r b + 3r 2 a rt + 15r 2 b r£ - 5rg) In s ab , mo( 

(63r^ + 90r 3 r ab - 62r 2 r 2 b - 30r a r 3 b - r 4 ab + 90r 3 r b - 90r 2 r ab r b + 30r o r 2 b r b 
b r b - 126r 2 r b - 90r a r ab r 2 - 62r 2 b r b + 90r a rj) + 90r ab r 3 + 63r b ) 

+ ^ (—3^a + 2r 2 r 2 b + ri b + §r 2 a rl + 2 r 2 ab r 2 b - 3r£) In s ab , 

= 7 Q 55 Q 0 ( 531r a + 6 30 r l r ab - 531r^r 2 b - 420r 3 r 3 b - 531r 2 r^ b + 630r a r 3 b + 531r® b + 630r 3 r b 

- 630 rir ab r b - 630 r a r 4 ab r b + 630r 3 b r b - 531r"r 2 - 914r 2 r 2 b r 2 - 531r" b r 2 - 420r 3 r 3 - 420r 3 b r 3 

- 531r 2 r b - 630r a r ab r b - 531r 2 b r b + 630r a rf + 630r ab r^ + 531r b ) 


A L r a r b := 


A 2 r a r b := 


3600 
~30r 3 b r b 
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3360 


- 3r® + 3r 4 r 2 b + 3r 2 r 4 b — 3r® b + 


3r a r b 




A 


-3 


Vb := 


1219276800 
- 31500r a r 7 b - 35825r® b 


(l8975r® + 18900rlr ab - 21100r®r 2 b - 18900r 3 r 3 b - 50550r 4 r 4 b + 65100r 3 r 3 b + 88500r^r« b 

6300r 3 r 2 b r b + 6300r 4 r 3 b r b - 44100r 3 r 4 b r b 


- 18900r 7 r b — 18900r®r ab r b 

- 4 - 2 - 2 ■ 7560r 3 r 3 b r 2 


- 44100r 2 r 3 b r b + 31500r a r® b r b - 31500rl b r b - 18204r®r b 2 + 3780r 3 r ab r b 2 - 30804r 4 r 2 b r b 2 

- 39492r 2 r 4 b r b + 44100r a r 3 b r 2 + 88500r 3 b r 2 - 8820r®rg - 3780r 4 r ab r 3 - 840r 3 r 2 b r 3 - 7560r 2 r 3 b r 3 

- 44100r a r 4 b r b + 65100r 3 b r b 3 - 1542r 4 r£ - 3780r 3 r ab r 4 - 30804r 2 r 2 b r^ + 6300r a r 3 b r 4 - 50550r 4 b r£ 

- 8820r 3 r b 5 + 3780r 2 r ab r b 5 - 6300r a r 2 b r b 5 - 18900r 3 b r b 5 - 18204r 2 r® - 18900r a r ab r b 6 - 21100r 2 b r® 

+ 18900r a r b + 18900r ab r b 7 + 18975r b 8 ) + EE( - 15r 8 + 20r 3 r 2 b + 30r 4 r 4 b - 60r 2 r 3 b + 25r® b + 12r 3 r 2 

+ ^rlrltrl + 36r 2 r 4 b r b - 60r® b rg + 6r 4 ^ + 12r 2 a r 2 ab r% + 30r 4 b r£ + 12r 2 rg + 20r 2 b r£ - 15rf) lns ab , 

(B28i) 


r 3 

A -1 — 

n 


Y^g( - 63r a - 90r 3 r ab - 2r 2 r 2 b - 90r a r 3 b ■ 

-I- 126r 2 r b + 90r a r ab r b + 126r 2 b r b - 90r a r 3 - 90r ab r b 
1 

+ 40 


63r 4 b + 150r 3 r b + 90r 2 r ab r b + 90r a r 2 b r b + 150r 3 b r b 


- 63r b ) 


( 3r a + 2 r 2 r 2 b + 3r 4 b — 6 r 2 r 2 — 6 r 2 b r b + 3 r b ) lns ab . 


(B28j) 


Appendix C: Explicit results for the field functions 

In this appendix we give the explicit formulas for the 
PN approximate solutions of both the constraint (ESI) 
and the field (14.31) equations which can be found in the 
literature or were computed by us for the first time. The 
different inverse Laplacians needed to compute field func¬ 
tions presented below are given in Appendix [Bl 

1. Potentials 0( n ) 

The solutions of Eqs. (13.101) for the potentials 0( 2 ) 
and 0( 4 ) are known in d dimensions. They can be ob¬ 
tained by means of Eqs. (lB2l) and (lB9l) using the result 
0 (2))reg(x a ) = m b r 2 ~ d . They read 

0 ( 2 ) = K^2m a r 2 a ~ d , (Cla) 


0(4) = 

, d — 2 

(Clb) 

2* (4)1 + 4(d~ lp 4)2 ’ 

5 (4)1 = 

^ m a 

a 

(Clc) 

S(4)l = 

a b^a 

(Cld) 


where the coefficient n is defined in Eq. (18.41) . These 
solutions are valid also for general n-body, i.e., not only 
for two-body, point-mass systems. 

The potential 0(g) is split into two parts, 0( 6 ) = 
0 ( 6)1 + 0 ( 6 ) 2 , where 0(6)i is the solution of the Poisson 
equation (13.14al) and 0 ( 6)2 fulfills the equation (17.71) . The 
solutions to these equations are fully known only for two- 
body point-mass systems and in d = 3 dimensions (they 
can be found in an implicit form e.g. in Refs. |82j,|83( and 
f75l]). They can be symbolically written as 


0(e) 1 " A 1 (E (-^4) + mJa 



(C2a) 




1 

87T 


1 



(C2b) 
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A 


more explicit form of the function 4>[q)\ can be written as@ 

1 __ y^ (Pg) 2 _ 1 y- \" WbP Q ( A I J_ 

32t t ^ mlr a (16tt) 2 ^ rn a r ab \r a n 

_ 9 1 y^ 2p% — (n a • p a ) 2 1 y^ y' I 

" 4(1 67 t) 2 4- r a (16-) 2 V&1 


+ 

1 

4 


_J_ yyrrftm, / j_ + J_\ 

(!H g r afe V r “ 

(Pa • V a ) (p b ■ V 6 ) (V a • V b ) 2 (A _ 1 r a r b ) 


+ [-8 (Po • Pfe) (Va • v b ) + 3 (p a • V a ) (p 6 • V b ) - 8 (p a • V b ) (pb • V a )] ( A 1 — ] 

V Vb ) 

+4 (Pa ■ Va) (Pb • Va) (V Q • Vft) ^A _1 | , 


(C3) 


where the inverse Laplacians A 1 (r a rb), A 1 (r a 1 r b x ), and A 1 {r a r b 1 ) can be found in Eqs. (1B28I) . The Poisson 
integral needed to calculate the function 4>{&)2 reads (here b ^ a) 


A- 1 



1 fflqWb 
32 (16tt) 2 rlnrlh 



12 ry ab 


l %r 2 a rl b 


- 12 r a r 


3 

ab 


4- 3 r 


4 

ab 


+28 rln - 12 r 2 a r ab r b - 12 r a r 2 ab r b + 28r^ b r b - 30r^r^ + 60r a r a br 2 - 30 r 2 ab rl - 36 r a r b - 36 r ab rl + 35r^) 


3 Pa 8 (n Q • Pa) 2 
647r m a r ^ 


1 1 

327T 77Tb 


/ 2 (p 6 - Vb) 2 - 
l ’’a 



r a r b 


2Pb (V • Vb) 2 - 4 (p 5 • V) ( P b • Vb) (V • V 6 ) 


-4 ( P b • V) ( P b • Vb) (V • Vb)] (a- 1 -) + 8 ( P b • V) 2 (V 1 —) + \ ( P b • Vb) 2 (V • Vb) 2 ( &- 1 -) } , (C4) 

\ Ta J \ Ta'T'b / U \ ^a/ J 

where the inverse Laplacians A _1 (rbr“ 1 ), A^ 1 (r“ 1 tV 1 ), and A -1 (r? r" 1 ) are given in Eqs. (1B28I) . 


2. Longitudinal field momenta 7r)U 

Equations (|3.16D fulfilled by the longitudinal field mo¬ 
menta 7 can be written in the form 

€)j= s m’ ^ 

where S^ stands for the source term. Making use of the 
decomposition (EHD one rewrites Eq. (1C 5 II in terms of the 
vectorial function VA^, 

A ^) + ^r d a v U = S U- ^ 

It is not difficult to find the formal solution of Eq. m 
in terms of the first and the second inverse Laplacian of 
the source S). ,. It reads 

+, = A- 4A W >s|, r (C7) 


5 In Appendix E] we employ the following notation: V = ( di ), 
V a = {dai), (Pa * V) = PaA, (Pa • Vfo) = Pai^ti, (V a * Vfe) = 
dai&bi • 


Let us now assume that the source term in Eq. m lias 
the form of a divergence of a symmetric and trace-free 
quantity, 

SU = 9^, (C8) 

where = Tj\ and T ^ = 0. Then, making use of 
Eq. (12.91) and the definition (12.141) of the TT operator in 
d dimensions, one can show that 

^) = t h-( t w) TT ( C9 ) 

Using the source term from Eq. (I3.16al) we have ob¬ 
tained the explicit solution (IC7l) for the leading-order 
function VA in d dimensions, 

= g(d-l) ^ ( (3d “ 2)Pai 

+ {d— 2) 2 (n a • pa) r4jr 2 ~ d . (CIO) 

This solution is also valid for n-body point-mass systems. 
The leading-order longitudinal field momentum can 

be constructed from the function V+ by means of the 
formula (12.91) . Let us also quote the following useful form 
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of the field momentum fr^ in d = 3 dimensions, 


= Ito 




^ ^ Pafc'j 2 
1 


&ik [ ] 4 “ &jk ( 

r a ) tj \ r a 


1 


a /,k 


2 r a,ijk 


(Cll) 


The next-to-leading-order function VA has been calcu¬ 
lated in d = 3 dimensions and for two-point-mass systems 
only. It reads 


= - 


16(167t) 2 ^ r\ 

1 


E + 2(n a • p Q )<) 


(167r) 2 


cl b^a 


EE m b < 4<9 ai (p Q • V) 


- 2 d i (Pa ■ v a ) + 4 p ai (V ■ V a ) 


A 


-1 


r a n 


- 2di (p a ■ V) (V ■ V Q ) A 


.-2 


1 


- d ai ( Pa • V a ) (V ■ V a ) [ A 


r a n 

r, 


n 


+ \di ( Pa ■ V Q ) (V ■ V Q ) 2 ( A" 2 ^ 
4 \ n 


(C12) 


where the inverse Laplacians A 1 (r a 1 r b 1 ), 


A 2 ( r a lr fc 1 ) I A 1 ( r a,r b 1 ), and A 2 {r a r b 1 ) are 
given in Eqs. (IB28I) . The field momentum n’A can be 
constructed from the vectorial function by means of 
the formula (EH). 


3. and related terms 


The leading-order TT part hj^ of the metric is the 
solution of Eq. (14.111) . It is the sum 


I TT _ J TT , J TT 
' t (4 )ij ~ ''(4)0 ij T" "(4)2 iji 


(C13) 


where hj^ 0i j does not depend on particles’ momenta 
and is quadratic in momenta. The quadratic-in¬ 

momenta part of hj^ 2i j can be computed in d dimensions 
and for general n-body systems. It reads 

h V)2ij = 8 ( ^_^ E^{[ p «-( rf + 2 )( n - • Pa) 2 ]Sij 
+ [d(d- 2)(n a ■ p a ) 2 - (d + 2)p 2 ]nX 


+ 2d(n a • p a )(n l a paj + n J a p a i ) + 2p ai p aj }i 


, 2 -d 


(C14) 

The momentum independent part hJ^ oi j is known only 
in d = 3 dimensions and for two-body systems. Its fully 
explicit form reads (here s ab := r a + r b + r ab for a ^ b; 
see, e.g., Ref. 


. TT 

n (4)0ij 


^EE«<! - ” ^< b < b + 2 + £)»X 


8 (167t) 
+ 16 


a b^a 
2 1 

S a& r lh 


Sab \^ab Sab 


' ab 


3 ab 


r a / r a 


r ab \ r b 


-1 - 


( n >a6 + KKb) 
17 4 


1 


r ab r a \r a 
1 4 


n . 8/1 1 

+ 3r a )-I-1- 

Sab \^a S a i ? 


« 


r ab r a r a r b s ab \r a r ab 


(C15) 


The next-to-leading-order conservative TT part of the metric function hj? is very complicated. It fulfills 

equation (14.121) . The piece of which diverges linearly at infinity in d = 3 dimensions equals A^ 1 h■ The part 
of this inverse Laplacian which corresponds to the function hj^ij can be computed in d dimensions (and for general 
n-point-mass systems). It reads 


a-Rtt _ 

' V ( A\ On n - 


d 2 


„4-d 


E —— { (16 - d 2 ){ n Q • p a ) 2 - (14 - d)p 2 S xj + 2(7 d - 8 )p a 

TYin v 


aiPaj 


d (4pij 48(4 _ d ^ d _ j) g t 2 Z_, ma 

+ (4 ~d) (4 + d)p 2 - (d - 2) 2 (n a • p a ) 2 « - 2(4 - d)(2d - l)(n Q • P a){n l a p a j + «>«*)} . (C16) 


The part of the inverse Laplacian A d corresponding to the function ' 0ij . can be computed only in 3 dimensions 
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and for two-point-mass systems. It reads 


/\-i /? TT 

A ft (4)0 ij 


1 d 2 x - v- f 1 

> , > , m a m b 

a b^a 


(167r) 2 dt 2 


1/3 2 \ 

-r - r a r 6 ) + — 

.'aft 


Sij + (<9, ; <9 b: ,- + djd bi - did a j - djd ai ) — 

^ ab 


+ uk\ didj ^ ~ r “ (r * + 17r “ fc) ] “ ^ ( A 1 ^) + ^ 


- W ( A~ 2 


1 


r a r& 


• (C17) 


The inverse Laplacians A _1 (r“ 1 rf -1 ) and A _2 (r“ 1 r^ 1 ) can be found in Eqs. (IB28I) . 

Other parts of the function hj^ are determined by A^ 1 Cj^ - and they enter the densit y h 2 '^ from Eq. (18.81) . We 
list below formulas by which, together with the inverse Laplacians enumerated in Appendix IB 21 one can compute the 
density h^ 2 ) explicitly in 3 dimensions. These formulas read 


W(3)) TT = 

4m 2 , 

(16tt) 2 1 

TT \ TT 

4 to 2 to 2 

'(4)ij ) 

(167t) 3 


( - 2{pudij +pijdu)^— + 

\ nr 2 2 


dudij (pi • Vi) — j + (l -o- 2), 


(C18a) 


(^)A^) - ^0 + 2 ^fS 13dlidlj ^ +15d2id2j o) + f 1 dudij w~ 2 


TT 


+ 


— (hpld 2 id 2 j - 4 (p2j&2i +P2id 2 j) (P2 • V 2 ) 'j —— + d 2 id 2 j (p 2 • V 2 ) 2 — 
(167r) z m2\^ / r\V2 T\ 


TT 


+ A-™( p|fcA^i 2 ') TT + (1 „ 2), 


127T m 2 

XT to 2 to 2 / 3 
L>/ K \„„ = 


n 2 


(C18b) 


3 (6)*7 (167t) 3 nr 2 

1 


1 35 


o dudij o 


TT 


4 rn 2 2 (~ ~ 1 


PI dijd 2i 


nr 2 


TT 


r\r 2 J (167r) 2 77j-i 

+ ^6^)2 (( 8 (Pi • P2 )dud 2j + IQpu^dij (p 2 • V 2 ) - d 2j (p 2 • Vi)) - 8 pup 2j (Vi ■ v 2))^- 
+ (d u (pi ■ Vr) (d 2j (p 2 ■ Vi) - pu (p 2 • V 2 )) + p 2j d u (pi ■ Vi) (V x • V 2 )) ^ 

+ (pu (P 2 • v 2 ) ^3 d 2 j (Vi • V 2 ) — -9yA 2 ^ — d\id 2 j (pi • V 2 ) (p 2 ■ V 2 )^ — 

+ -^u (pi • Vi) (p 2 • V 2 ) (~d 1 jA 2 - d 2 j (Vi • V 2 ))rir 2 j 

x TT 

d-, + h j + (1 -H- 2). (ci 8 c) 

1 - 


r 1 2 x TT 

5 m 2 pupij 1 PiPiiPy 

- 0 1 -o 

327T m 1 ri 2 4 m\ 


To compute the density explicitly one has to employ 
Eqs. (IC18I) together with the following inverse Lapla¬ 
cians, which can be found in Eqs. (iBbfl) and (IB281) : 


Let us also note that after combining and 


(6 )ij 

V TT 


A-M a , 

A~ 2 5 a , A ~ 3 S a , 

A -1 1 

a — 2 ^ a —: 

rir 2 

.La , x_A 

nr 2 

r 2 

A' 2 — , A -3 —, 

r 2 r 2 

A _1 nr 2 , 

A~ 2 nr 2 , A -3 : 


nr 2 


(^( 2 )A/iTT.) tt i nto (B (6)ij - + (1/4) <f> (2 )Ah^ ..) 11 , see 
Eq. (18.911 . there is no need to compute any in¬ 
verse Laplacians of (rfr 2 )~ 1 . The inverse Laplacian 
A -1 (<(>( 2 ),kihj^ kl ) also needed to compute h^ 12 ) explic¬ 
itly, can be inferred from Eqs. (IC2bl) and (IC4I) . 


4. Local d-dimensional UV analysis 

In Sec. 3 of Ref. [3l| a method of local analysis of UV 
divergences in d dimensions was proposed. This method 





































48 


was essentially used in [3lJ for studying the local behavior 
of only one specific function, namely, the momentum- 
independent part hJ^ 0i j of the field function At 

the 4PN level we have had to use it in many more cases; 
therefore, we present this method here in more detail. 

The problem is to find the local behavior in d dimen¬ 
sions and, say, around x = xi, of the solution of equation 
of the type (for some cases we have to consider the Pois¬ 
son equation without TT projection) 

(A d f) TT = g TT - (C20) 

We assume that the full analytic solution of this equation 
is not available. The source function g depends on the 
field point x only through x — Xi and x — X 2 ; therefore, 
one defines the auxiliary function g as follows: 

g{n) := 3 (x-xi,x-x 2 ) = g(nni,nni +ri 2 ni 2 ), 

( c 2i) 

where we have used x — xi = nni and x — x 2 = 
rini + ri 2 ni 2 . The method of Ref. [3l| relies on expan¬ 
sion of the source function g around rq = 0 and applying 
the operator A^ 1 to each term of the expansion. This 
expansion has the form 


g(n) = 


V g (fc) ( n i) k 

^ k\ 15 


(C22) 


where to > 0 is some nonnegative integer. With the 
expansion (IC22D of the source term one can relate the 
following expansion of the solution to Eq. (1C20I) near 
x = xi: 

OO 

fl mhom( X ) = A rf 1 (3 (fe) ( n l) r l) TT - (C23) 


Let us write the formal solution of Eq. (IC20I) in the 
form of the integral 

/ TT ( X ) = -« J d d o; , g(x , )(|x — x'| 2 ~ d ) TT . (C24) 


The crucial element of the method is expansion of the 
kernel of the integral (IC24I) around x = xi. To do this 
one introduces the auxiliary function 

K(n) = |x - xf ’~ d = |(x - Xl ) - (x' - Xl )| 2 - d 

= |nni - r' 1 n[f~ d , (C25) 


and expands it around rq = 0: 


K(n) = 


nj.rj) 


’=0 


l\ 


(C26) 


One then substitutes the expansion (IC26I) into the inte¬ 
gral (IC24I) and integrates the sum term by term. Conse¬ 
quently one gets the expansion of the form 


°° 1 r 

/hol( x ) = / d d x’ g(x’)(K w (n i;n' 1 ,rj)rf) TT . 

(C27) 

As an example of applying the above formula let us 
compute the leading-order term in the expansion of the 
momentum-independent part of the function 

near x = Xi. Using Eq. (IC15I) one easily checks that the 
function is finite for x = Xi in d = 3 dimensions. 

This finite value cannot be obtained in d dimensions from 
the expansion (IC23I) , but it follows from Eq. (IC27I) (with 
i = 0). After making use of Eq. (I3.20bl) [which is the 
source term for the function , see Eq. (14.1111 ] one 
gets 


The above formula does not contain all terms of the 
expansion of the physically acceptable solution of Eq. 
(fC2fll) . The missing terms are solutions of homogeneous 
Poisson equation. We have devised a method to compute 
these terms. 


lTT 
Z (4)0 ij 


(x = Xl ) = 


(d + 1 )d(d — 2) 3 4 ac 2 
16(d- 1) 3 (4 -d) 


■TOi?U 2 


x ( 5ij — dn\ 2 «i 2 ) r 


.4-2 d 


(C28) 
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